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Abstract 
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In this dissertation, we aim to predict the evolution of the various physical properties of the binary 
semiconductor SrS after its doping with iron (Fe) and sp-type metals. In this context of study, we 
consider ternary compounds (SrS mono-doped with iron) and quaternary compounds (SrS co-doped 
with iron and sp-type alkali metals), some of which are studied in volume (3D) and others in monolayer 
form (2D). Properties calculations are based on the Full Potential Linearized Augmented Plane Wave 
(FP-LAPW) method and on the PBE and PBE + mBJ approximations implemented in the WIEN2K 
simulation program. The work involves assessing the impact of elements added to the SrS 


semiconductor host matrix on the properties of final compounds produced by doping and co-doping. 


Initially, we investigated the structural, mechanical, electronic and magnetic properties of single- 
doped Sri.xFexS compounds in the rock-salt structure, considering the following Fe concentrations: 0%, 
12.5%, 25%, 50% and 75%. The study is carried out at 3D and the main objective is to investigate new 
potential half-metallic ferromagnets (HMF) for spintronic applications. The results obtained show that 
the compounds containing 12.5%, 25% and 50% iron are half-metallic ferromagnets with a total 
magnetic moment value equal to 4u1, and are thermodynamically and mechanically stable; nevertheless, 


the iron-rich compound (75% Fe) is metallic. 


We then turned our attention to co-doping the SrS matrix to produce quaternary compounds. The 
latter are realized by doping the ternary compound (12.5% Fe) with alkali metals (Li, Na and K). The 
aim is to see the effect of the sp metals on the properties of the ternary compounds. The study of these 
new compounds, also in 3D, is based on density functional theory (DFT) and semi-classical Boltzmann 
theory (BT), and covers structural, electronic, magnetic, optical and thermoelectric properties. The 
results obtained show that the incorporation of alkali metals into the SrS: 12.5% Fe compound renders 
them half-semiconducting (HSC) with narrower energy gap values than those of HMEF ternary 
compound (SrS: Fe). The Curie temperature (T.) values obtained are higher than room temperature. In 
addition, the magnetic moment of quaternary compounds is higher than that of ternary compounds, 
having a value of 5ug. Examination of the optical properties revealed a shift of the spectra towards the 
visible region, accompanied by a broadening of the compounds' absorption band. Investigation of the 
thermoelectric properties of p-type compounds showed that their figure of merit (ZT) values are greater 
than unity at 1200K, underlining their exceptional transport efficiency and making them highly 


promising candidates for optoelectronics and high-temperature thermoelectric applications. 


The 2D study focused on the structural, electronic and magnetic properties of ternary compounds. 


The calculations are based on a plane-wave pseudopotential (PAW) method using the PBE and 


II 


PBE+HSEO06 functionals, integrated into the VASP simulation software. The effect of dimensionality 
reduction was mainly observed at the level of electronic structures, where the character of ternary 
compounds becomes HSC instead of HMF observed in their 3D counterparts. This result has a positive 
impact on the optical and transport properties of ternary compounds. Under these conditions, the total 


magnetic moment resulting from p-d hybridization is 4ug per unit cell. 


Keywords: Ab-inito calculations, SrS, Half-metallic ferromagnetic, Alkali co-doping, 3D-bulk 


materials, 2D-monolayers, Spintronics, Optical properties, Thermoelectric properties. 
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Résumé 


Dans cette thése, nous visons a prédire |’évolution des différentes propriétés physiques du semi- 
conducteur binaire SrS aprés son dopage au fer (Fe) et aux métaux de type sp. Dans ce contexte d’ étude, 
nous considérons des composés ternaires (SrS mono-dopé au fer) et des composés quaternaires (SrS co- 
dopé au fer et aux métaux alcalins de type sp) dont certains composés sont étudiés en volume (3D) et 
d’autres sous forme de monocouche (2D). Le calcul des propriétés est basé sur la méthode d'ondes planes 
augmentées linéarisées a potentiel plein (FP-LAPW) ainsi que sur les approximations PBE et PBE+mBJ 
implémentées dans le programme de simulation WIEN2K. Le travail consiste 4 évaluer l'impact des 
éléments ajoutés a la matrice semi-conductrice SrS sur les propriétés des composés finaux réalisés par 


dopage et par co-dopage. 


Dans un premier temps, nous nous sommes intéressés aux propriétés structurales, mécaniques, 
électroniques et magnétiques des composés mono-dopés de type Sri..FexS dans la structure rock-salt en 
considérant les concentrations de Fe suivantes : 0 %, 12,5 %, 25 %, 50 % et 75 %. L’ étude est faite a 
3D et l'objectif principal est d'investiguer de nouveaux ferromagnétiques semi-métalliques (HMF) 
potentiels pour des applications dans le domaine de la spintronique. Les résultats obtenus ont montré 
que les composés contenant 12,5%, 25% et 50% de fer sont demi-métalliques ferromagnétiques avec 
une valeur du moment magnétique total égale a 4ug et sont thermodynamiquement et mécaniquement 


stables ; néanmoins, le composé riche en fer (75% Fe) s’aveére métallique. 


Dans un deuxiéme temps, nous nous sommes intéressés au co-dopage de la matrice SrS pour réaliser 
des composés quaternaires. Ces derniers sont réalisés par le dopage du composé ternaire (12,5% Fe) par 
les métaux alcalins (Li, Na et K). Le but est de voir l’effet des métaux sp sur les propriétés des ternaires. 
L’étude de ces nouveaux composés, toujours a 3D, repose sur la méthode de la théorie de la fonctionnelle 
de la densité (DFT) et sur la théorie semi-classique de Boltzmann (BT) et concerne les propriétés 
structurales, électroniques, magnétiques, optiques et thermoélectriques. Les résultats obtenus montrent 
que l'incorporation des métaux alcalins dans le composé SrS : 12,5% Fe les rend demi-conducteur (HSC) 
avec des valeurs du gap énergétique plus étroites que celles du composé ternaires HMF (SrS : Fe). Les 
valeurs de la température de Curie (T.) obtenues sont supérieures a la température ambiante. En outre, 
la valeur du moment magnétique des composés quaternaires est supérieure a celle des composés 
ternaires, elle vaut 5g. L'examen des propriétés optiques a révélé un décalage des spectres vers la région 
du visible accompagné d’un élargissement de la bande d’absorption des composés. L’étude des 
propriétés thermoélectriques des composés de type p a montré que les valeurs de leurs figures de mérite 


(ZT) sont supérieures a l'unité a la température 1200K, soulignant l'efficacité exceptionnelle de leur 


IV 


transport et faisant d’eux des candidats trés prometteurs pour l'optoélectronique et pour les applications 


thermoélectriques 4 haute température. 


L’étude a 2D a concerné les propriétés structurales, électroniques et magnétiques des composés 
ternaires. Les calculs effectués reposent sur une méthode de pseudo-potentiel 4 ondes planes (PAW) en 
utilisant les fonctionnelles PBE et PBE+HSE06, intégrées dans le logiciel de simulation VASP. L'effet 
de la réduction de la dimensionnalité a été principalement observé au niveau des structures 
électroniques, oti le caractére des composés ternaires devient HSC au lieu de HMF observé dans leurs 
homologues a 3D. Ce résultat impacte positivement les propriétés optiques et de transport des composés 
ternaires. Dans ces conditions, le moment magnétique total résultant de I'hybridation p-d, vaut 4up par 


cellule unitaire. 


Mots-clés : Calculs ab-inito, SrS, Demi-métaux Ferromagnétiques, Co-dopage alcalin, Matériaux en 


bulk-3D, Monocouches-2D, Spintronique, Propriétés optiques, Propriétés thermoélectriques. 
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General Introduction 


he inception of the transistor in 1947 sparked a pivotal turning point in the evolution of 
communication and technology [1]. With the ability to digitize information on-chip, the world 
entered the information age, where knowledge and communication became the driving forces 
of economic prosperity, surpassing dependence on physical effort and natural resources. In today's 
world, we witness persons communicating through cell phones, working on laptops, and utilizing the 
World Wide Web for various tasks. Integrated circuits and high-performance electronics have become 
integral to our lives, powering applications in drones, recording devices, contemporary cars, and 
manufacturing machinery, and more. This remarkable growth and evolution of technology has 
revolutionized the global economy, profoundly affecting all aspects of everyday life, science, industry, 


and technology. 


The continuous growth of technology-based semiconductors has been driven by the ever-growing 
data processing capabilities, rapid integration density, and faster speeds. Advancements in magnetic 
materials and optical have further facilitated the seamless transmission, processing, and storage of vast 
amounts of information. Despite these impressive achievements, the semiconductor-based technology 
is now facing significant challenges. According to the International Technology Roadmap for 
Semiconductors (ITRS), further downsizing and power reduction will be hindered as the miniaturization 
of elements size and operational speed approach their limits [2]. Additionally, electronics inherently 
generate waste heat during switching, leading to higher power consumption in electronic devices. To 
overcome these obstacles, the focus is on developing innovative particle-less technologies to handle the 
growing information demands. Simultaneously, researchers are exploring new ways to efficiently 


manage and organize information in this rapidly evolving technological landscape. 


Exploring alternative integrated circuits presents a promising path forward, and one such example 
is the realm of all-optical devices utilizing photons as information carriers - the fundamental quanta of 
electromagnetic waves. This exciting field has given rise to photonic crystals, a subset of optical 
materials characterized by their periodically adjusted refractive index [3-5]. Just as a semiconductor's 
periodic atoms lattice creates energy bands and bandgaps for electrons, photonic crystals introduce 
frequency bands and bandgaps for electromagnetic waves, known as photonic bandgaps. When light (or 
photons) falls within these bandgaps, the photonic crystals efficiently reflect it, regardless of the incident 


angle. What makes photonic crystals truly fascinating is their ability to control and manipulate their 
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properties. By modulating physical parameters such as the dielectric constant or lattice parameters, the 
characteristics of the photonic bandgaps can be finely tuned. This unique selectivity has resulted in 
various applications, from photonic waveguides to integrated circuits, unleashing the potential for 


advanced optical devices [6-10]. 


In recent times, nano-magnetism and magnetization dynamics [11] have emerged as compelling 
subfields within the realm of spintronics [12], capturing significant attention by virtue of their promise 
in revolutionizing nanoscale signal processing and information exchange mechanisms [13-15], like 
filters, transistors, and storage elements [16-18]. Spintronics, originally abbreviated as "SPIN Transport 
Electronics," revolves around transmiting, processing, and storing informations hinged on the 
magnetization state of a given system. A pivotal breakthrough in spintronics occurred with the discovery 
of the giant magnetoresistance (GMR) effect in Fe/Cr multilayered systems by Fert and Grunberg [19, 
20], for which they were honored with the Nobel Prize in 2007. The GMR effect, dependent on the 
relative orientation of two ferromagnetic layers when subjected to a spin-polarized current, 
revolutionized magnetic hard disk drives and paved the way for novel non-volatile magnetic memories 
like Magnetoresistive Random Access Memory (MRAM) [13]. Since then, additional discoveries, such 
as Spin Transfer Torque (STT) [21], facilitating current-assisted magnetization switching, and the 
(Inverse) Spin Hall Effect ()SHE) [22], enabling the generation and detection of spin currents, have 
further propelled the field. The Spin Pumping Effect [23], the Spin Seebeck Effect (SSE) [24], and the 
Dzyaloshinskii-Moriya Interaction (DMD) [25] have also introduced fascinating possibilities to the world 


of spintronics. 


Besides spintronics, the SSE gives rise to the nascent scientific field of so-called spin-caloritronics, 
in which transport design methodologies have been considered as future independent power sources 
since they enable the generation of electricity from waste heat in magnetic materials, thus positioning 
it supplementary to the well established domains of spintronics and thermoelectricity [26]. The 
field of spin-caloritronics, derived from the term "calor" (Latin for heat), and which came into existence 
following the significant discovery of the spin Seebeck effect (SSE) in 2008 by Uchida, Saitoh, and their 
colleagues [24] emerges from the intricate coupling of spin, charge, entropy, and energy transport in 
predominantly magnetic structures and devices [27]. This quickly developing sector finds several uses 


in coolers, power generators, and thermometers [28]. 


In light of these compelling factors, the pursuit of advanced and efficient spintronic, photonic and 
thermoelectric materials gains even greater significance for the emergence of eco-friendly energy 
devices in the realm of materials science and technology. Diluted Magnetic Semiconductors (DMSs) 
that are semiconducting materials doped with magnetic impurities, typically transition metals (TM), and 


that demonstrate room temperature ferromagnetism, are widely employed in spintronics due to their 
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remarkable half-metallic ferromagnetic nature (HMF), exhibiting a 100% spin-polarization in both their 
bulk form (3D) and monolayer form (2D) [29]. Interest in the utilization of DMSs is further fueled by 
the prospect of finely tuning the bandgap through different processes such as doping, co-doping, 
alloying, etc., thereby providing the ability to precisely control light absorption at specific energy levels, 
tailored to suit the desired application [30]. Another important feature is the high carrier mobility 


resulting from the quantum confinement [30], which is widely desirable for thermoelectric devices. 


The extensive findings and in-depth comprehension of the physical characteristics of these materials 
owe much to ab-initio calculations utilizing Density Functional Theory (DFT). This computational 
approach minimizes the time, risks, and expenses associated with experimental studies. Furthermore, 
the results obtained through ab-initio calculations generally align with experimental observations. At 
the very least, such calculations can predict material properties, thus providing valuable guidance for 
experimental investigations. A notable example is in the field of spintronics, where the captivating 
property of half-metallicity, as documented in the literature [29], was initially revealed through ab-initio 
calculations [29], subsequently confirmed through experimental validation [31]. In the realm of photonic 
and optoelectronic materials, the groundbreaking concept of "valleytronics" in 2D materials, which 
enables precise control of light emission and absorption through manipulation of the valley degrees of 
freedom, has also taken the same path [32]. In recent times, ab-initio calculations have experienced 
remarkable advancements, primarily due to significant progress in computational techniques and 


computing power. 


Given the significance of DMSs materials and their potential applications in the various desired 
domains, we employ, in this dissertation, ab-initio calculations to systematically investigate spin- 
dependent structural, mechanical, electronic, magnetic, optical, and thermoelectric properties of selected 


host matrix semiconductor SrS-based DMSs using different processes. 


The layout of the dissertation is structured as follows: 


Chapter 1 entitled “Literature Background”, serves as concise introduction to the domains of 
spintronics and spin caloritronics, providing an overview of diluted magnetic semiconductors (DMSs) 
that includes a diverse array of materials, spanning from 3D bulk semiconductors to low-dimensional 
2D materials. Additionally, the chapter delves into a comprehensive bibliographic research on the 
fascinating realm of doped and co-doped structures, where elements ranging from alkali to transition 
metals are integrated into the host SrS to amplify the different physical properties of the materials under 


investigation. 
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Chapter 2 entitled “Theoretical and Computational Formulation”, is dedicated to an in-depth 
exploration of the Density Functional Theory (DFT), which forms a fundamental basis for the ab-initio 
calculations. It delves into the PAW (Projector Augmented Wave), LAPW (Linearized Augmented 
Plane Wave), and FP-LAPW (Full Potential Linearized Augmented Plane Wave) methods, elucidating 
their significance in our research. Various approximations used in our study are thoroughly examined 
and discussed. Moreover, this chapter provides a comprehensive overview of the functionalities of the 


WIEN2K and VASP codes, both of which were crucial tools employed in our research. 


Chapter 3 and 4 


In each of chapters 3 and 4, a comprehensive compilation of the calculation parameters, the obtained 


results, and their insightful interpretations is presented. 


Chapter 3 entitled “Effect of Doping and CO-Doping on the SrS Bulk Properties”, focuses on 
investigating the spin-resolved structural, mechanical, electronic, and magnetic properties of bulk rock- 
salt Sr).xFe,S at various concentrations, namely x = 0, 0.125, 0.25, 0.50, and 0.75. Additionally, the 
chapter delves into the structural, electronic, mechanical, magnetic, optical, and thermoelectric 
properties of bulk SrS: Fe co-doped with alkali metals Li, Na, and K, while maintaining a fixed 
concentration of x = 0.125. The chapter further compares the properties of these co-doped materials to 


those of singly doped Sro.g7sFeo.125S. 


Chapter 4 entitled “Effect of Dimensionality Reduction on Fe-Doped SrS Properties”, deals with the 
structural, electronic, and magnetic properties of monolayers SrjxFexS at various concentrations, 
including x =0, 0.125, 0.25, 0.50, 0.75, and 1. The chapter further compares the properties of these 


doped-monolayers to their bulk counterparts. 


Finally, the dissertation concludes with a comprehensive summary covering the most important 


findings and future perspectives. 
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Chapter 1 


Literature Background 


his chapter provides a concise yet comprehensive introduction to the fields of spintronics and 
spin caloritronics as well as dilute magnetic semiconductors (DMS), which include various 
3D bulk semiconductor groups as well as low-dimensional 2D materials. Furthermore, it 
explores the fascinating area of doped/codoped structures, in which alkali to transition metals are 
incorporated to improve structural, mechanical, electronic, magnetic, optical and thermoelectric 
properties. By exploring the fascinating properties and diverse uses of these materials, this chapter helps 


in selecting suitable compounds for this dissertation. 
1.1 Introduction 


The ever-increasing energy demand of human activities has become an urgent global problem and 
requires the introduction of advanced techniques for research and development of sustainable and 
environmentally friendly energy resources. The scientific and industrial community is extremely curious 
to find an ideal renewable energy source that can be practically translated into the development of useful 
products and mitigate the harmful effects of anthropogenic carbon emissions on the global climate. 
Among the various options, the thermoelectric phenomenon stands out as an excellent source of 
renewable energy. As electronic devices continue to shrink toward the nanometer and their operating 
speed increases, dissipating heat or waste energy has become a critical issue. To achieve energy savings 
and improve the performance and reliability of electronic devices, it is imperative to either reduce or 
effectively utilize wasted energy. In this context, spintronics proves to be a promising path to 
significantly lower energy consumption [1]. On the other hand, thermoelectrics explores the direct 
conversion of waste heat into electrical power and represents a fascinating avenue for research and 
development. A fascinating and novel research area known as spin caloritronics is currently gaining 
significant attention as it combines the strengths of spintronics and thermoelectrics [2]. Semiconductor 
materials, whether in their bulk form or as 2D counterparts, with their intrinsic properties and well- 
suited electronic structures emerge as perfect candidate bases for this purpose. Based on their energy 
bandgap, semiconductor materials are divided into three different categories, of which we are interested 


and will only discuss the wide bandgap semiconductor (WBG) category. 
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1.2 Spintronics and Related Devices 


1.2.1 Definition and Historical Development of Spintronic Phenomena 


Spintronics (or spin electronics) has revolutionized electronic data storage by actively manipulating 
the spin degrees of freedom of electrons interacting with their orbital moments, ushering in a new era of 
possibilities [3]. Magnetic layers act as spin polarizers or analyzers and determine spin polarization 
through spin-orbit coupling. This complicated interaction leads to the generation of a spin current, which 


is further driven by fascinating spin waves. 


Every electron possesses two possible states known as spin-up <> and spin-down & , Meaning 
they can rotate either clockwise or counterclockwise. These majority-up and majority-down domains 
are randomly distributed. However, when an external magnetic field is applied, it aligns the domains in 
the direction of the electric field. This phenomenon enables precise control and manipulation of electron 
spins and forms the basis for various fast and energy-efficient spintronic applications and innovations 


in electronic data storage and processing [4]. 


The attractive field of spintronics began in 1922 with the discovery of the magnetic moment of 
electrons [5]. Wolfgang Pauli's groundbreaking research validated the quantization of electron spin and 
introduced the idea of Pauli matrices. A crucial moment came in 1973 when the conversion of electron 
spin to spin current between ferromagnetic films was observed [6] and shortly afterwards D'yakonov 


and Perel predicted the spin Hall effect [7], which was later confirmed experimentally. 


In 1975, Julliere and his team conducted the first experiment on the tunneling magnetoresistance 
effect (TMR), revealing the fascinating phenomenon of electron tunneling between two ferromagnets 


separated by a thin insulator [8]. 


In 1976, the concept of generating spin-polarized current in a semiconductor by passing a current 


through a ferromagnet/semiconductor junction was first proposed [9]. 


In the 1980s, the idea of spintronics emerged, which explores spin-dependent electron transport in 
solid-state devices. The journey continued with the discovery of spin-polarized electron injection from 
one ferromagnetic metal to another by Johnson and Silsbee [10]. This was followed by the independent 
giant magnetoresistance (GMR) breakthrough, observed in ferro/metal/ferro multilayer structures, 
reported by Fert et al. and Griinberg et al. Their pioneering work earned them the Nobel Prize in Physics 
in 2007 [11, 12]. 


A milestone in the practical application of spintronics was the introduction in 1997 of spin valve 


sensors in read heads of hard disk drives, which were developed by an IBM researcher [13]. This 
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breakthrough enabled an impressive storage capacity of 16.8 GB, setting the stage for other 


manufacturers to follow suit and continually improve performance and storage capacities. 


However, despite remarkable progress, spintronics has faced several challenges. The search remains 
for low-cost, abundant, and lightweight ferromagnetic materials with 100% spin-polarized current, 
longer spin lifetime, high magnetic anisotropy, and easy manipulation of spin currents. Furthermore, it 
is important to consider the influence of strain and temperature on spin properties to exploit the full 


potential of spintronic devices [14]. 


During the late 20" century, heavy transition metals such as iron (Fe), nickel (Ni), and cobalt (Co) 
gained recognition as important components of spintronic materials due to their high Curie temperatures 
(TC) and intrinsic ferromagnetic (FM) nature. These properties make them ideal for generating and 
controlling spin polarization, which is fundamental to various spintronic applications. Due to their 
remarkable contributions, these materials have laid the foundation for remarkable advances in this field 


[15]. 


As a result, the field of spintronics experienced rapid growth, driven by the growing interest and 
contributions of numerous research groups. This trend can be seen in the significant increase in 


publication numbers over the years, as shown in Figure 1. 


Number of Publications from Dimensions Database 
Search Key : Spintronics 


2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 


-e- Publications (total) 


Figure 1. Annual publications on spintronics over the past decade. 
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1.2.2 Spin Devices 


Figure 2 provides a comprehensive and insightful overview of spintronic devices, which can be 
conveniently divided into two distinct categories: Mott-type and Dirac-type, with each category 
characterized by its unique use of electron spins, spin waves, and spin/orbit moments. In the Mott-type 
devices, the most important phenomena are giant magnetoresistance (GMR) and _ tunneling 
magnetoresistance (TMR), while the Dirac-type devices exploit the fundamental spin-orbit interactions 


[16]. 


To better understand the evolutionary progress of these spintronic devices, we can divide them into 
three generations, each with a notable leap in their capabilities and applications. The first generation 
revolves around spin transport, with electrical spin generation playing a central role. The second 
generation includes devices that exploit spin dynamics, making skilful use of spin-orbit effects, electric 
fields and electromagnetic wave applications. These advances have significantly expanded the horizons 
of spintronics and led to transformative functionalities and innovative possibilities. Finally, the 
groundbreaking developments of the third generation pushed the boundaries of spintronics even further 
to include three-dimensional structures and quantum technology. This development has catapulted 
spintronics into the realm of quantum computation and promises a revolutionary paradigm shift in the 


field. 
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Figure 2. Catalog of spintronics devices. After Ref [16]. 
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1.2.3 Spin Generation 


The generation of spin-polarized electrons in non-magnetic (NM) materials can be achieved by 
various methods [16]. These techniques include spin injection from ferromagnetic (FM) materials, 
application of magnetic or electric fields, introduction of electromagnetic waves, Zeeman splitting, spin 
driving force, thermal gradients, and mechanical rotation (see Figure 3). A commonly used approach is 
spin injection through FM materials such as conventional FM metals (Fe, Co, and Ni), half-metallic 


ferromagnets (HMF), and dilute magnetic semiconductors (DMS). 


Besides, spin-polarized electrons in an NM material may exhibit a population difference due to a 
stray field near the edge of an FM. Electromagnetic waves such as circularly polarized light and 
microwaves can also excite spin-polarized electrons in SC materials, where the optical selection rule 
plays a crucial role. Conversely, a spin-polarized electron current can produce circularly polarized light 
emission. These principles can be extended to spin generation by other electromagnetic waves, including 


spin pumps and high-frequency spin induction. 


Furthermore, the application of a thermal gradient due to the spin-Seebeck and Nernst effects has 
been shown to be effective in generating spin-polarized carrier currents and offers valuable opportunities 


for energy harvesting and other potential applications. 


Electrical spin 
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thermal gradient wave 
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Figure 3. Methods for generating spin-polarized electrons in non-magnetic media. After Ref [16]. 


11 


Chapter 1 Literature Background 


1.3 Opportunities at the Frontiers of Spintronics 


As mentioned in the previous section, an effective method for generating spin-polarized electrons 
is through the application of a thermal gradient, which leads to the emergence of a new field at the 


forefront of spintronics, known as “spin caloritronics”. 


Spin caloritronics deals with the study of thermoelectric transport effects in material systems in 
which magnetic spin moments are introduced. The inclusion of additional spin degrees of freedom not 
only opens up a wealth of conceptually innovative mechanisms and functionalities, but also ushers in a 
new era of transformative breakthroughs in the conversion of thermal to electrical energy in solids. In 


Figure 4 we present a clear representation of the concept of spin caloritronics. 


Spin 


Caloritronics Spintronics 


(ep 


Thermoelectrics 


Figure 4. An illustrative depiction of the concept of spin caloritronics. 


Since its founding in 2009 through an international workshop, the field of spin caloritronics has 
attracted considerable attention. Annual meetings of specialists in the field have become the norm as the 
research topic attracts more and more researchers. In particular, the twelfth “Spin Caloritronics” 


workshop took place in 2023 [17]. 


The German Physical Society (DFG) has recognized the importance of this emerging field and 
initiated a priority program to support further research [18]. In addition, the US Department of Energy 
has provided funding for the SHINES (Spins and Heat in Nanoscale Electronic Systems) research 


cluster, which includes 14 research groups from 7 different institutions across the country. 


The interdisciplinary nature of spin caloritronics, combining magnetism, thermoelectrics and 
microelectronics, ensures that it will continue to captivate the global scientific community and generate 


immense interest in cutting-edge research [19]. 
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Before continuing our discussions in this new area, it is first necessary to look at the history of 


thermoelectricity to better understand the role of spin in our study. 


1.3.1 Thermoelectric Effects 


Thermoelectricity is an intrinsic property of materials that gives them the ability to convert thermal 
energy into electrical energy. Over time, this fascinating area has experienced two significant periods of 


development: 


The first stretched from 1821 to 1851, the second lasted from the late 1930s to the early 1960s. In 
1821, the pioneering German physicist Thomas Johann Seebeck made a significant discovery, revealing 
one of the fundamental thermoelectric effects [20]. His observation involved a metal needle placed 
between two different metal conductors connected at their ends by junctions at different temperatures, 
resulting in a noticeable deflection. He originally incorrectly attributed this phenomenon to a magnetic 
field induced by the temperature differences at the metal junctions and suggested that this could explain 


the earth's magnetic field. 


However, in 1825, Oersted corrected this misunderstanding by showing that the effect was due to 
the emergence of a potential difference at the junction of two materials with temperature contrast. This 
discovery became known as the “Seebeck effect.” Consequently, this effect allows the use of a 


temperature gradient to generate an electric current, as described by the equation: 


AV = Sap. AT (eq.1) 


With Sp, being the Seebeck coefficient or thermoelectric power of the two materials A and B, typically 
expressed in uV.K”, and AV representing the electric potential difference generated by the temperature 


difference AT. 


In 1834, physicist and watchmaker Jean-Charles Peltier made the groundbreaking discovery that a 
temperature difference occurs at the junctions of two different materials when they are exposed to an 
electric current [21]. However, it took until 1838 for Lenz to fully explain the phenomenon [22]. 
Through a convincing experiment, Lenz showed that the direction of heat transfer depends on the 
direction of current flow. He crystallized water around a bismuth-antimony compound and by reversing 


the direction of the current; he managed to melt the ice. This phenomenon became known as the “Peltier 


effect.” 
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Since each material has its own Peltier coefficient zr, flowing a current through the junction results 
in an interruption in the flow of heat at that point. For example, as current flows from material a to 


material b, the junction heats up (if a > 7b) or cools (if a < mb), as described by the equation: 
Q = Myp.1 (eq.2) 


With Q representing the generated thermal power, I denoting the electric current, and T,4p being the 


Peltier coefficients of the two materials. 


In 1840, a crucial breakthrough occurred when James Prescott Joule revealed his groundbreaking 
discovery. He showed that whenever an electric current flows through a material, an amount of heat is 
generated that is directly proportional to the intensity of that current [23]. This remarkable phenomenon 


became known as the “Joule effect” and its expression is captured by the following equation: 


qj = p.J? (eq.3) 
With p representing the electrical resistivity of the material. 


It is obvious that the square of the current flux density and the positive nature of the electrical 
resistivity contribute to the inevitable positivity of the amount of heat, q;. Unlike the reversible Peltier 
and Seebeck effect, the Joule effect is irreversible as it merely produces the production of heat without 
the ability to absorb it. In addition, the physicist William Thomson established the crucial connection 
between the Seebeck and Peltier effects in 1851 [24]. Consequently, when a temperature difference and 
an electric current are applied at the same time, a dynamic heat exchange with the environment occurs. 
Conversely, a material exposed to a thermal gradient and with heat flowing through it produces an 


electric current, a process known as the “Thomson effect.” 


The heat output or absorbed per unit volume can be expressed using the following formula: 


ot 


Q; = -t.J.VT (eq.4) 


Tt denotes the Thomson coefficient. The Thomson coefficient is closely linked to the Seebeck and Peltier 


coefficients, as expressed in the following relationship: 


TAR = Sap-T (eq.5) 
ds 
t= Pp. a (eq.6) 
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The Peltier, Seebeck and Thomson effects have a fundamental difference: the latter phenomenon 
manifest itself within a single material, eliminating the need for a connection between different 


materials, as in the case with other two. 


The timeline of discovery of these various thermoelectric effects by their respective pioneers is 


shown in Figure 5(a), along with schematic representations of each effect (Figure 5(b)). 
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Figure 5. (a) Key pioneers in thermoelectricity. (b) The representation of Seebeck, Peltier, and 


Thomson effects, respectively. 
1.3.2 Transport Equations 


Thermoelectric devices are based on the connection of two pairs of materials, one of type p with S 


> O and the other of type n with S < 0. When electric current is applied, charge carriers move from the 
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cold to the hot source. More specifically, electrons in the n-type branch migrate toward the hot source, 
while holes in the p-type branch also migrate toward the same heat reservoir. This concerted action 
results in both charge carriers transporting entropy from the cold source to the hot source, ultimately 


inducing a heat flow that actively counteracts heat conduction [25]. 


Within each branch, the total flow is expressed by: 


aT 


Qp = SpIT — KpAp = (eq.7) 
dT 
Qn = SIT — Ky An (eq.8) 
Where Kp and K, are the thermal conductivities of the p-type and n-type materials, Ap and Ay are their 


respective cross-sectional areas, Sp) and S, are the Seebeck coefficients, and z represents the spatial 


coordinate. 


Heat is dynamically transported from the cold source to the hot source, resulting in a total flux Q; 


[25]: 
Q:= (Q, os Qn)|,-o (eq.9) 


At the same time, the Joule effect comes into play in the circuit due to the flow of the electrical 
current I. The heat generated by this effect is given by pl?/A. Conservation of energy can be aptly 
described for both branches of the circuit, taking into account the balance achieved by a non-constant 


thermal gradient: 


ar Ip 

KyAp dat = re (eq.10) 
aT fp, 

An Gr = rik (eq.11) 


When analyzing the system, we must carefully consider the boundary conditions. Let Ln and Lp 
denote the lengths of the individual branches. At the cold source (z = 0), the temperature remains the 
same as the cold source, ensuring a consistent thermal interface. Similarly, at the hot source (z = Lp or 


z= Ln), the temperature is maintained at the hot source level, thereby maintaining thermal equilibrium. 


The boundary conditions can therefore be formulated precisely as follows [25]: 


T=T, at z=0 (eq.12) 
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T=T, at z=L, orz=L, (eq.13) 
The eq.10 and eq.11 yield: 
aT I pp(2—5Lp)  kpAy(Tr-Te) 
KpAp ie = a ar + a a (eq.14) 
aT Ppp (2-5Ln) Ky An(Th-T¢) 
KnAn a = = a a ae (eq.15) 


By substituting these formula into eq.7 and eq.8 and employing eq.9, we reveal the comprehensive 


expression for the total thermal flux Q:: 
1 
Qt = (Sp — Sn) IT ¢ — KAT — 5°R (eq.16) 


Where K is the thermal conductivity and R is the electrical resistance of the circuit, both defined as 


follows: K 


_ KpAp | KnAn 

k= i. + is (eq.17) 
_ LpPp , LbnPn 

R= “Ay + ae (eq.18) 


The combined Joule and Seebeck effectsplay a crucial role in determining the dissipated power W 


[26]: 
W =1.[(Sp — Sy). AT + IR] (eq.19) 


The efficiency & of the thermoelectric cooler is defined as the ratio between the extracted heat Q; 
and the dissipated electrical power W. By utilizing equations eq.16 and eq.19, it follows [26]: 
Q: _ (Sp-Sn)IT¢—KAT-H?R 


Eo = 


WT [(Sp—Sp)-AT+IR] (eq.20) 


In addition, it is possible to determine the efficiency of a p-n device specifically designed for 
electricity generation from a temperature difference. The efficiency &c is defined as the ratio of the useful 


electrical power delivered to a load resistance r and the thermal flux through the device: 
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= Wuseful = I[(Sp —S;,).AT+IR] 


Ee Q 172 
(Sp Sy )IT s—KAT-5I (R+r) 


(eq.21) 


This profound expression provides a way to discover the electrical current that maximizes 
efficiency. Both cooling and power generation have two different values of electrical current I, each 


optimizing either the conversion efficiency or the electrical energy generated by using heat. 


When both efficiencies reach their peak, an interesting observation occurs that shows their 
dependence exclusively on the temperatures T. and T; and the dimensionless figure of merit Zpn.Ta 
[27]. Here Ta = (T. + Ty)/2 represents the average temperature, a crucial parameter in this thermoelectric 
context. The figure of merit factor Zpn.Ta is a critical factor that is defined specifically for the particular 
pair of materials used in the device. The Zpn.T formulation includes the intrinsic absolute parameters 


of the materials constituting the thermoelectric couple [27]: 


(Sn-Sp)* 


(a ee a) 
C/ KpPpty KnPn)” 


T (eq.22) 


pn 
The maximum occurs when the efficiency reaches its peak. 


Similarly, individual p-type and n-type materials each have their own intrinsic factor, commonly 


known as the figure of merit: 


ZT = —T =—T= T (eq.23) 


The figure of merit factor ZT is of utmost importance when assessing the thermoelectric properties 
of a material and serves as a crucial indicator of its suitability for thermoelectric applications. This 
dimensionless factor plays a critical role in whether a material has desirable thermoelectric properties 
or does not meet the criteria. In order to achieve the greatest possible efficiency of thermoelectric 
conversion, a maximum value of ZT must be aimed [28]. This requires maximizing the Seebeck 
coefficient and achieving the highest possible electrical conductivity and hence the numerator S7o, 
called the power factor, while at the same time minimizing the denominator, which represents the sum 


of the electronic contribution to thermal conductivity, ke , and the lattice contribution, «. 


1.3.3 Optimization Parameters for ZT 


We have just seen the importance of the quality factor in identifying high-performance 
thermoelectric materials, whether they belong to the n-type or p-type semiconductors. To realize the 
greatest potential of ZT, these materials must have distinct and precise thermal and electrical transport 


properties. However, we face a challenge because these transport phenomena have strong correlations. 
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Notably, a material with excellent thermal conductivity often also has exceptional electrical conductivity 


[29]. 


The insightful Figure 6 illustrates the overall evolution of thermoelectric properties with respect to 
catrier concentration and provides valuable insights into the complex relationship between thermal and 


electrical properties in thermoelectric materials. 


Figure 6. Evolving thermoelectric parameters with charge carrier concentration, n, at 300 K. After Ref 


[30]. 


This figure provides a compelling insight into the complex interdependence of thermoelectric 
parameters on carrier concentration. As the concentration of charge carriers varies, a delicate 
equilibrium is created in which the Seebeck coefficient decreases while the electrical and thermal 
conductivity increases. It becomes clear that there is an optimal carrier concentration that maximizes the 
figure of merit typically found in heavily doped semiconductors. Remarkably, this optimal concentration 


is in the range of 10!” to 10° cm®, although each semiconductor has its own optimum. 


When exploring the fascinating material properties described in the literature, a promising 
thermoelectric material should manifest itself as a glass due to its low thermal conductivity to phonons 
and as acrystal to electrons due to its exceptional electrical properties [31]. Achieving peak performance 
requires a complex crystal structure that hinders the propagation of phonons. Cutting-edge research has 
identified numerous degrees of freedom, allowing researchers to influence the figure of merit and shape 


material performance [32]. These critical factors are illustrated in Figure 7. 
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Figure 7. Degrees of freedom affecting the figure of merit. After Ref [30]. 


These four degrees of freedom make it possible to optimize the quality factor. These mechanisms, 
which are crucial for charge carrier transport, are the key to understanding the complex interactions 
between charge, orbitals, spin, and lattice structure [33]. For example, by manipulating the crystal 
complexity of the semiconductor, the proportion of optical phonons increases compared to acoustic 
phonons. Since acoustic phonons mainly transport heat, a more complex structure leads to reduced 
thermal conductivity [34]. The incorporation of impurities or point defects increases phonon scattering 
and effectively reduces thermal conductivity [35]. At the same time, doping improves the electrical 
properties of the semiconductor [36]. Furthermore, heat transport through dimensionality reduction to 
limit the phonon mean free path and approximate the interatomic distance follows a mode of random 
thermal diffusion known as Einstein's model [37]. Furthermore, dimensionality reduction leads to 
significant fluctuations in electron density, favoring a high Seebeck coefficient without affecting the 


electrical conductivity of the material [38]. 


1.3.4 Boltzmann Transport Theory 


The simulation of electrical conductivity and the Seebeck coefficient can be achieved by applying 
first principles calculations alongside the use of the Boltzmann transport equation. Charge transport 
occurs in the presence of an electric field and/or a temperature gradient. This phenomenon is described 


in the following equation [39]: 


J=ed f.v=o0E (eq.24) 
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In this context, J represents the flow of charge; e stands for the electronic charge; f denotes the charge 


distribution; o denotes the electrical conductivity; E represents the electric field, and v corresponds to 
the charge velocity. By understanding the temporal and spatial distribution of charges, we can 


determine the flow of charges. The temporal development of the charge distribution can be illustrated 


as follows [40]: 
d 
Fy e Vf + 2Vpf = es 2 (eq.25) 


Here, r represents the electron's position; p signifies momentum, and c denotes collision index. More 


specifically, eq.25 illustrates the change in charge distribution after a collision. 


fe ia) 
oo (eq.26) 
at 
f-fo=Cerz Tt (eq.27) 
By utilizing eq.25, eq.26, and eq.27, we can derive: 
0 
f=fote (- a TV.E (eq.28) 


Upon substituting eq.24 with the expression provided in eq.28, the conductivity can be effectively 


described as follows: 
) 
o =e (yu? (eq.29) 
This equation can be reconfigured into tensor form as part of electronic structure calculations: 


o(€) = = =f de — (2) es kTn, kUn, kUn, KO(E — En k) (eq.30) 


In this context, Q signifies the volume of the unit cell; e denotes the charge of the carrier; € represents 
the energy of the group; N stands as the count of k-points employed in the computation; fg corresponds 
to the Fermi-Dirac distribution function; t denotes the relaxation time, which represents the average 
time elapsed between two collisions; v signifies the charge group velocity, and 6 is the delta function 


[41]. The indices k and n refer to the crystal momentum and band index, respectively. 


The velocity,uv, can be estimated from the band structure through the following relationship: 
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(eq.31) 


Here, # is the reduced Planck constant. 


In the presence of a temperature gradient, the total electric field deviates from the expression in 
eq.24 due to the influence of the Seebeck effect. This change results in the adjusted formulation of eq. 


24 as shown below: 
J = oE — oSVT (eq.32) 
Moreover, the heat flux (Q) generated from the temperature difference can be defined as follows: 
Q=TJ, = ST] — KVT (eq.33) 
Where Js represents the entropic flow. Eq. 32 and eq. 33 are recognized as Onsager relationships [42]. 


Using the Onsager equations, we can derive the Seebeck coefficient and electronic thermal 


conductivity from band structure calculations: 
_ “2 on} Ofo 
S= Jde a Ge Ge me) nk Te kUn, kUn, KO(E — Enk) (eq.34) 


2 
kpT 4 


Ke = WO oO if de — (2) (= ay Sn kTn, kUn, kUn KO(E —€Ey k)- Tos? (eq.35) 


Here, w denotes the chemical potential, and kg stands for the Boltzmann constant. 


Based on the above equations, we can define electrical conductivity, Seebeck coefficient and 
electronic thermal conductivity. The two unknown parameters are the lattice thermal conductivity and 


the value of the relaxation time. How to determine them is already discussed in Chapter 3. 


1.3.5 Spin-Powered Thermoelectrics 


With thermo-spin effects, energy transfer is controlled by the spin current instead of the 
conventional charging current. This fascinating phenomenon has revealed several novel thermo-spin 
conversion processes in magnetic materials and their complex interconnecting structures. Notably, 
notable phenomena such as the spin-dependent Seebeck-Peltier effect [43, 44] have come to the fore, 
offering enormous potential for revolutionizing spintronics and applications in thermal energy 


conversion. 
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The spin-dependent Seebeck-Peltier effect is a fascinating thermo-spin conversion phenomenon in 
ferromagnetic materials (FMMs). Its essence arises from the inequality of the Seebeck-Peltier 
coefficients between electrons with spin-up and spin-down orientation. To understand the transport 
phenomena within FMMs, a two-current model that treats spin-up and spin-down electrons as different 
entities is required [45]. The electron density of states (DOS) within FMMs depends on spin, leading to 
spin-sensitive variations in transport properties, including thermopower, electrical conductivity, and 


electronic thermal conductivity. 


Using the two-current model, a comprehensive picture of the transport properties is created, which 
includes the entire electrical conductivity, the Seebeck coefficient and the electronic part of the thermal 
conductivity. This is achieved by combining the relevant information from the two spin channels, as 


given in the following equations [45]: 


For electrical conductivity, the formula is as follows: 


o=ott+oal (eq.36) 


Accounting for the complicated nature of the spin-Seebeck coefficient requires careful 
consideration of the up- and down-spin coefficients in relation to their respective electrical 


conductivities: 


S= (oT ST +oalS1)/(oT +01) (eq.37) 


The electronic part of thermal conductivity can be expressed by the following simple summation: 


ke = keT +Kel (eq.38) 


However, the lattice thermal component works independently of spin considerations. We use a 


separate model to calculate this aspect. Further details will be explained in a later discussion. 


In today's applications, the Peltier and Seebeck effects are of central importance for power 
generation and cooling. Semiconductors are playing a pioneering role in exploiting these phenomena. 
At absolute zero temperature (O K), these materials exhibit electrical properties similar to those of 
insulators. However, as temperature increases, electrons acquire significant ability to generate electric 
currents [46]. Essentially, semiconductors occupy an intermediate position, having electrical 
conductivity levels that span the spectrum between metals and insulators. A semiconductor is classified 
as intrinsic if it maintains a pure state, although achieving absolute purity is virtually unattainable. In 
such cases, charge carriers arise exclusively from crystal defects and thermal excitation [29]. 


Consequently, these semiconductors have minimal electrical conductivity and have limited application 
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in thermoelectric scenarios. To address this challenge, doping becomes a crucial technique by 
introducing impurities into the intrinsic semiconductor matrix. In addition, the semi-metallic property, 


which is vital for spintronics, can only be realized through the process of doping semiconductors. 


Subsequently, our goal will shift to the classification of wide bandgap semiconductors, a special 


category that is the focus of our research and piques our great interest. 
1.4 Wide Bandgap Semiconductors : Host Base for Different Applications 


Wide-gap (WBG) semiconductors are a fascinating class of semiconductor materials characterized 
by their relatively large energy band gap (E, > 2 eV). Due to their exceptional properties, these materials 
have become essential components of today's electronic devices and energy applications. They feature 
adjustable electrical conductivity, controlled carrier concentration and high optical transparency, 
making them extremely versatile and desirable. Among the prominent wide bandgap semiconductors, 
transparent conductive oxides (TCOs) have been intensively studied since the 1950s [47]. Materials 
such as tin-doped indium oxide (ITO) and amorphous In-Ga-Zn-O (IGZO) are extensively used in 
display technologies and solar cells [48]. These TCOs exhibit a remarkable combination of transparency 
and conductivity, allowing their use as transparent electrodes, enabling efficient light transmission while 


enabling effective electrical conduction. 


While much of the research and device integration efforts have been directed toward wide bandgap 
oxide materials, it is important to note that the scope of transparent semiconductors extends beyond 
oxides. Since the early 20" century, numerous classes of non-oxide semiconductors have been 
experimentally demonstrated to exhibit both transparency and conductivity. In this context, the 


following classes have acquired particular importance: 


* Carbides with group IV anions [49] such as silicon carbide (SiC) and nitrides with group V anions 
[50] such as gallium nitride (GaN) have revolutionized power electronics due to their excellent 
electronic properties. SiC and GaN offer significant advantages, including high breakdown voltages, 
excellent thermal stability, and impressive electron mobilities. These properties make them ideal for 
power conversion and control applications, enabling more efficient energy management and improved 


performance in a range of electronic systems. 


* Halides like y-Cul with anions from group VII [51] and IV two-dimensional (2D) materials [52] 
like graphene have shown great promise in various optoelectronic devices. These materials have unique 
properties that can be exploited for applications such as advanced sensors, ultra-fast electronics, and 
next-generation photonic devices. Their exploration opens up exciting possibilities for the development 


of innovative and powerful optoelectronic technologies. 
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* Furthermore, extensive research efforts have been dedicated to group-based VI chalcogenide (Ch = S, 
Se, Te) semiconductors, encompassing various material compositions. This includes binary II-VI MCh 
semiconductors, such as ZnS [53], and CdS [54], as well as other binary M,Chy semiconductors (e.g., 
SnSz2, In2S3, where M represents a metal). Ternary chalcopyrite I-II1]-Chz semiconductors, predominantly 
based -copper, have also received significant attention. Examples of this class include CuAIS>2 [55], as 
well as other ternary compounds like a-BaCu2S2 [56] and Cu3TaS4 [57]. Moreover, multinary-layered 
mixed-anion compounds, such as LaCuOCh [58], have been investigated for their unique properties. 
Additionally, research has extended to 2D chalcogenides, both binary and ternary materials, including 


MoS:, further broadening the scope of exploration in this field. 


The particular chemistries and properties of WBG chalcogenide semiconductors distinguish them 
from other compounds, particularly their oxide counterparts, and offer specific advantages for 
optoelectronic applications and as p-type semiconductors. This distinction can be understood by 
considering factors such as the covalency of metal-VI (M—VI) bonds, which tend to increase as we 
move down group VI in the periodic table. This increase in covalency leads to larger orbital overlaps 
and lower hole effective masses. Heavier chalcogenides, with elements such as sulfur (S), selenium (Se), 
and tellurium (Te), possess higher-lying p orbitals (S-3p, Se-4p, Te-5p) that can hybridize with metal 
(M) orbitals, such as copper (Cu) 3d orbitals. This hybridization results in the formation of more disperse 
and delocalized valence bands compared to oxides [59]. Consequently, p-type chalcogenides exhibit the 
potential for achieving higher mobilities than their oxide counterparts. For example, Cu-based 
chalcogenides have demonstrated p-type mobilities of up to 20 cm? V! s! [56], while Cu-based wide 
bandgap oxides typically exhibit mobilities lower than 1 cm? V's". 

Figure 8 provides a comprehensive network diagram architecture highlighting various WBG 
material classes, with particular emphasis on chalcogenides, which are clearly shown in a brown 


gradient. The classes are systematically organized and sorted according to their respective anion groups. 
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Figure 8. A network topology diagram illustrating the diverse material categories within WBG 


semiconductors. After Ref [60]. 


On the other hand, the simplest class of WBG chalcogenide semiconductors in terms of their 


respective cation groups consists of binary divalent metal chalcogenides designated M*?Ch” (Ch = S, 


Se, Te). Among these, the II-VI chalcogenide materials (referred to as II-Ch) are the most commonly 


used binary materials for electronic and optical applications. They typically contain cations of groups 


IIA and IIB, including elements such as Mg, Ca, Sr, Ba, Zn or Cd, among others. 


The family of strontium chalcogenides SrX (X = S, Se, Te) has recently attracted considerable 


attention in both experimental and theoretical studies due to their intriguing physical properties [61, 62, 


63]. Their chemical composition leads to a number of fascinating properties, including the ability to tune 


their band gaps, exhibit high electrical conductivity, exhibit exceptional thermal stability, and display 


excellent optical properties spanning a wide spectrum of electromagnetic waves. In addition, strontium 


chalcogenides offer significant advantages compared to other materials, making them extremely 


attractive for scientific research and technological applications. Their abundance in nature, cost- 


effectiveness and environmental friendliness contribute to the growing interest in harnessing their 


potential for various advances in science and technology. 
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Furthermore, the versatile SrX chalcogenide family offers a compelling opportunity for extrinsic 
doping with both isostructural and non-isostructural binary systems. This opens up a wealth of 
possibilities that enable the exchange of cations and anions as well as the formation of ternary and/or 
quaternary compounds within the common crystal structure. A particular focus in research on SrX 
chalcogenides is the addition of extrinsic elements to the matrix in different concentrations [64—68]. 
This has become a fascinating area of research with the aim of modifying the inherent indirect nature of 
the electronic band gap and improving the optical absorption properties. By reducing the large band gap, 
these materials can be individually tailored to the specifications of optoelectronic devices. This trend 
was inspired by the successes observed in related technologies such as CIGS (copper indium gallium 
selenide) and CdTe (cadmium telluride), where the integration of extrinsic species to passivate intrinsic 
defects, the introduction of flat dopants, and the use of substitution doping (alloying) for band gap 
technology has led to impressive solar cell efficiencies of over 22% [69]. The integration of extrinsic 
species into the SrX chalcogenide family also offers exciting prospects for state-of-the-art spintronic 


devices and efficient spin caloritronic systems [70-72]. 


From revolutionizing and advancing various fields to enabling efficient energy storage and 
facilitating the right path to energy conservation, research into strontium-based chalcogenides promises 


to address critical challenges and promote innovation in multiple areas. 


For this dissertation we chose strontium monochalcogenide (SrS). Recognizing the importance of a 
comprehensive review, the following section delves into an in-depth study of both pure and SrS-based 
dilute magnetic semiconductors (DMS). In particular, we investigate their properties and applications in 
three-dimensional (3D) bulk structures as well as two-dimensional (2D) monolayer structures to 


highlight their importance in the fields of spintronics, optoelectronics, and thermoelectronics. 
1.4.1 Bulk SrS Monochalcogenide 


Strontium sulfide (SrS), a promising member of the alkaline earth sulfide (AES) family, has proven 
to be a highly sought-after phosphor in various fields. Extensive research efforts have been devoted to 
exploring its applications in various fields, revealing its extraordinary potential. SrS has demonstrated 
superior performance and versatility in both bulk form and thin films, making it an attractive material 
for numerous practical applications. One of the main strengths of SrS lies in its luminescence properties, 
which have been extensively studied and exploited. Its indirect wide bandgap of 4.32 eV in bulk form 
[73] makes it an ideal choice for luminescence applications. This band gap allows SrS to create suitable 
luminescence centers and emit visible light without the undesirable effect of self-absorption [74]. This 
property ensures efficient light emission and improves the overall performance of SrS as a phosphor. 
The remarkable versatility of SrS phosphor is evident in its wide range of applications. In the field of 


optoelectronics, SrS is used in the development of efficient blue-emitting diodes (LEDs), displays and 
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radiation dosimetry. In addition, SrS is used in medical imaging to improve image quality and precision 


and enable more accurate diagnoses and treatments. 


In addition, SrS has applications in various other areas, such as lighting technology, where it 
contributes to energy-efficient and long-lasting lighting solutions. SrS's versatility extends to quantum 
technologies and advanced sensing systems, where its luminescence properties play a critical role in the 
development of next-generation devices. SrS exhibits polymorphism, meaning it can crystallize in 
different phases depending on temperature and pressure conditions. At room temperature and ambient 
pressure, SrS adopts a rock salt crystal structure, also known as B1 structure (Figure 9(a)). Alternatively, 
it can be represented as a regular face-centered cubic (FCC) lattice arrangement of anions, with cations 
occupying the octahedral holes in the lattice. Both types of ions (strontium cations Sr?* and sulfide 
anions S*) have a coordination number of 6 with an overall compactness of 74%. This makes it an open 
structure and offers the possibility of introducing both light and heavy atoms. Figure 9(b) also presents 
a comprehensive representation of the Brillouin zone belonging to the FCC Bravais lattice. In this figure, 


the irreducible part of the Brillouin zone is highlighted and surrounded by prominent orange lines. 


Figure 9. a) The Unit-cell of the semiconductor SrS, which adopts a rock-salt structure. Green spheres 
represent the Sr atoms, while orange spheres depict the S atoms. The image was generated using VESTA 
software [76]. b) The first brillouin zone of an FCC crystal, with the irreducible wedge enclosed by 
orange lines. The special points within the brillouin zone are identified by their common names, and the 


reciprocal base vectors are labeled as b, bz, and b3. After Ref [77]. 


In the B1 structure, each strontium cation is surrounded by six sulfide anions and vice versa, 
resulting in six-fold coordination. The FCC lattice ensures a symmetrical and densely packed structure 
and gives the crystal stability. This arrangement of atoms enables efficient charge transport and 


favorable electronic properties and contributes to the usefulness of SrS as a semiconductor material. The 
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six-fold coordination in the rock salt structure ensures the optimal balance between electrostatic 
interactions and the packing efficiency of the atoms within the crystal lattice, contributing to optimized 
electrical conductivity. This coordination arrangement influences the mechanical, optical, electronic and 
thermoelectric properties of the material and plays a crucial role in determining its behavior in various 


applications. 


However, fascinating phenomena occur in SrS under high pressures. Syassen et al. conducted X- 
ray diffraction experiments and discovered a pressure-induced structural phase transition in SrS at 
around 18 GPa [73]. This transition involves a transformation from the NaCl-type structure (B1) with 


six-fold coordination to the CsCl-type structure (B2) with eight-fold coordination. 


Furthermore, it was predicted that after the B1-B2 phase transition, the non-metallic nature of SrS 
would give way to metallic behavior at even higher pressures. This transition from a non-metallic to a 


metallic state is attributed to a mechanism known as bandgap overlap [75]. 
1.4.2 Tuning the Properties of SrS Monochalcogenide 


The properties of bulk SrS can be fine-tuned and optimized using a variety of methods, allowing 
properties to be tailored to specific applications. Doping [64-68], co-doping [78, 79], alloying [80], 
defect engineering [81], dimensionality reduction [82, 83] and external stimuli [84] are powerful 


techniques for manipulating and improving the properties of these binary connection. 


¢ Doping: involves the introduction of foreign atoms into the SrS lattice to change its properties. By 
selectively adding dopants, we can change the electronic and optical properties of SrS. Doping can 
influence the band gap, carrier concentration and conductivity, enabling the optimization of SrS for 


applications in optoelectronics, photovoltaics, thermoelectrics and sensing. 


* Co-doping: refers to the simultaneous introduction of two or more dopants into SrS. This technique 


allows for more precise control of material properties by combining the effects of multiple dopants. 


¢ Alloying: Solid solutions are created by alloying SrS with other materials, which expands the range of 
properties that can be achieved. By incorporating different elements, researchers can change the lattice 


structure, band gap and physical properties of SrS. 


¢ Defect engineering: involves the targeted manipulation of defects within the SrS crystal lattice. By 
controlling the concentration and type of defects, researchers can fine-tune the electrical, optical and 
mechanical properties of SrS. Defect engineering plays a crucial role in improving the thermoelectric 


efficiency, luminescence properties and mechanical strength of SrS. 
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¢ Dimensional reduction: refers to the reduction of bulk dimensionality of SrS into low dimensional 
structures such as 2D monolayers and their derivatives, i.e. nanotubes (1D), nanoribbons (1D) and 
nanoflakes (OD). At low dimension, SrS exhibits unique properties due to quantum confinement effects 
and an increased surface area to volume ratio. Reducing the dimensionality of SrS provides improved 
electronic and optical properties, improved charge transport, and increased surface reactivity, making it 


valuable for applications in nanoelectronics, catalysis, and energy storage. 


¢ External stimuli: such as temperature, pressure and electric fields can be used to modify the 
properties of SrS. By exposing SrS to certain environmental conditions, its physical properties can be 


changed. External stimuli also enable exploration of SrS behavior under different conditions. 


Although a lot of work has been done using these methods, especially in the area of doping and 
alloying, there is still a need for a systematic study of this compound and its derivatives using different 
approaches. In this dissertation, we chose doping, co-doping and dimensionality reduction as the primary 
techniques to tune the properties of rock salt SrS. These techniques have proven effective in improving 
material performance in optoelectronics, spintronics, spin caloritronics and nanoelectronics, which are 
indeed the future technological trends. The basic building blocks of spintronics, dilute magnetic 
semiconductors (DMS), were selected by doping with transition metals and studied in both 3D bulk 
configurations and 2D monolayers. Alkaline metals were also selected as co-dopants due to their three- 


dimensional volume structure. 


Figure 10 provides a representation of the extrinsic elements examined in this dissertation. The 
binary host elements are marked in green, the dopant elements in red, and the elements that form 
substitutional co-dopants are shown in blue. Certain elements such as Sr, S and Fe are considered for 


both 3D and 2D doping scenarios (in yellow). 
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Figure 10. Visualizing the covered elements in this dissertation within the periodic table. 


Moving forward, we will continue our discussion about DMS and alkali-substituted elements, 


highlighting their specific and exceptional properties as well as their applications. 


1.4.2.1 Diluted Magnetic Semiconductors (DMSs) 


Magnetic semiconductors are unique materials with both intrinsic semiconductor properties and 
magnetic properties. These substances provide a rich platform for theoretical and experimental 
investigations, allowing researchers to explore core ideas in condensed matter physics. The focus of 
magnetic semiconductor research is the complex interplay between semiconductor band electrons, such 
as charge carriers or impurity states, and the magnetic moments of the constituent ions. This unique 
interaction leads to a variety of electronic phenomena, including giant magnetoresistance, robust 
magneto-optical rotation, magnetically driven semiconductor-metal transitions, and fascinating 
magnetic polaron effects, among others [85]. These phenomena are missing from conventional 


semiconductors, making magnetic semiconductors a fascinating subject of research. 


While rare earth chalcogenide compounds, particularly those containing europium, have attracted 
considerable attention, the initial fascination with these materials stemmed from the discovery of EuO, 
a structurally simple semiconductor ferromagnet. At that time, EuO was only the second known 


semiconductor ferromagnet [86]. The ferromagnetic behavior of these compounds was explained by the 
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indirect coupling between local moments on the europium ions and the conduction band states. 
Nevertheless, the reported Curie temperature (Tc) of about 170 K, which is the temperature at which the 
metal loses its spontaneous magnetization - that is, the temperature at which the ferromagnetic phase 
changes to the paramagnetic phase - remains well below room temperature when it is endowed with Gd. 


This represents a significant obstacle in the search for improvements in this regard. 


Subsequent studies have focused their attention primarily on dilute magnetic semiconductors 
(DMS), which are currently among the most extensively studied systems. These DMSs are composed 
of alloys in which transition metal ions, including Cr**, Mn7*, Fe**, Co*, etc., are randomly incorporated 


in place of specific group II elements in II-VI semiconductor compounds, as explained in Figure 11. 
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Figure 11. Shematic illustration of different types of semiconductors based on their magnetic 
properties: (a) Magnetic semiconductor, (b) Diluted magnetic semiconductor, and (c) Non-magnetic 


semiconductor. Green circles represent the magnetic ion. After Ref [87]. 


On the semiconductor side, the multicomponent nature of DMS allows continuous tuning of 
structural and electronic parameters by adjusting the composition. On the magnetic side, they exhibit 
phenomena such as spin-glass transitions, antiferromagnetic short-range ordering, and anisotropic 
exchange interactions, making them fascinating random magnetic systems. Beyond their scientific 
appeal, DMSs hold promise for various technological applications. The ability to grow DMS in 
multilayer form with atomic precision using methods such as molecular beam epitaxy (MBE) opens 
opportunities for the integration of these materials into II-VI based integrated semiconductor structures 
[88]. This not only promotes their scientific research, but also enables practical advances in areas such 


as optoelectronics, spintronics, magnetoelectronics and spin caloritronics. 


The main goal of doping conventional semiconductors with magnetic ions is to make them 
ferromagnetic at room temperature without changing their semiconductor properties, and thus make 
them usable in spintronics. Although almost all materials have some degree of magnetism, only three 
are capable of maintaining permanent magnetization with Tc above 600 K: Fe, Ni and Co. These 


materials have high magnetic susceptibility, about 200 for Fe, and are called ferromagnetic. On a 
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microscopic scale, they consist of tiny regions in which each atom exhibits the characteristics of 
magnetic dipoles that are naturally oriented in a uniform direction. When exposed to an external 
magnetic field, these areas align with the field, resulting in significant overall magnetization of the 
material. Furthermore, this magnetization persists even when the external excitation is removed, a 


phenomenon known as hysteresis [89]. 


In fact, these materials represent typical ferromagnetic metals, but with a unique twist as they 
manifest as new semi-metallic ferromagnetic materials (HMF). The prediction of such materials dates 
back to 1983, when de Groot et al. [90] conducted a theoretical study on the band structure of the half- 
Heusler alloy NiMnSb. This alloy exhibited an intriguing property at the Fermi level (Er), where it had 
a zero electronic density of states for minority spins and a nonzero density of states for majority spins. 
Consequently, it conducts electricity for one spin direction (up), while for the other spin direction 
(down), it behaves like a semiconductor or insulator. This fascinating behavior allows achieving a 


remarkable 100% spin polarization, as indicated by the following relationship: 


tp agh 
— NertNer 
= ataciait 
Ner—Ner 


P (eq.39) 


Where: N io and N ee are the values of electronic density of states at the Er for spin-up and spin-down 


electrons, respectively. 


The difference in density of states between a non-magnetic metal, magnetic metal, non-magnetic 


semiconductor, and half-metal is schematically depicted in Figure 12. 
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Figure 12. Schematic representation of the difference in density of states between a non-magnetic 


metal, magnetic metal, non-magnetic semiconductor, and half-metal. After Ref [91]. 
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In certain DMS, interactions between electrons tend to align adjacent atomic magnetic dipoles 
antiparallel, resulting in two sublattices with opposite magnetization. The French physicist Louis Néel 
developed the first relevant model of this “anti-ferromagnetism” in 1936 [92]. In a “ferrimagnetic” 
material, neighboring atoms, due to their different natures, form magnetic dipoles with opposite 
orientations but different amplitudes (e.g. an Fe** ion and an Fe** ion), resulting in incomplete 
compensation. The spontaneous magnetization of such materials is non-zero below their Curie 


temperature [92]. Another visualization can be found in Figure 13. 
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Figure 13. Alignment of the magnetic spin of each individual iron atom. 


It is generally accepted that the presence of partially filled d or fbands is a prerequisite for a material 
to exhibit magnetism. These partially filled bands lead to complex exchange interactions that are crucial 
for the emergence of magnetic order. The exchange interaction, a cornerstone of quantum mechanics, 
plays a crucial role in determining the long-range magnetic order in ferromagnets [93]. However, its 
influence extends beyond ferromagnets and profoundly affects ferrimagnets and antiferromagnets, 
where neighboring magnetic ions are subjected to forces that align their individual moments into parallel 
or antiparallel configurations. This profound interaction results from the delicate interplay of Coulomb 
forces and the Pauli exclusion principle. Their effects can be effectively captured and are given by the 


Hamiltonian of the Heisenberg exchange as follows: 
Hey = —2) dizi SiS; (eq.40) 


In this context, S$; represents the spin operator of the i atom, and the parameter J is commonly referred 


to as the exchange integral. The above equation may be expressed as follows in the continuum model, 
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which takes into account the exchange interaction only between nearest neighbors inside the cubic 


lattice: 
Eox= Acx {(Vm)7dV (eq.41) 


Where, m = M/Ms is the continuous vector of reduced magnetization. The parameter Ag,, referred to as 
the exchange stiffness constant for the cubic lattice of spins, assuming the following form: 
_ 2ys? 


Aex = — (eq.42) 


a 


Where a is the lattice constant. 


The previously mentioned interaction is referred to as direct exchange interaction, which 
encompasses the interaction of electrons between magnetic atoms and their closest neighbors. However, 
this exchange can also take place indirectly, linking magnetic moments across relatively greater 
distances. There are several noteworthy types of indirect exchange mechanisms, such as the Ruderman— 
Kittel-Kasuya—Yosida (RKKY) exchange, which involves the coupling of metallic ions through itinerant 
electrons; super-exchange, occurring through various nonmagnetic ions mediating the exchange; and 
anisotropic exchange interaction (also known as Dzyaloshinskii-Moriya interaction) [94], where the 


spin-orbit interaction exerts significant influence and often results in a slight canting of neighboring 


spins. Additionally, there is the double exchange that arises between ions in different oxidation states. 


Other important interactions are the s-d and p-d exchange interactions, which refer to the interaction 
between the conduction and the valence bands, respectively, with the d electron states of the magnetic 
atom in a material. This interaction is crucial in confirming and reinforcing the ferromagnetic order in 


DMSs. 


The s-d and p-d exchange interactions in DMSs, arises from intentional doping of magnetic 
impurities (TM ions) into a non-magnetic semiconductor matrix. The s and p electrons of the host atoms 


interact with the d electrons of the magnetic impurities, creating local magnetic moments in the material. 
The Heisenberg exchange Hamiltonian in this case can be written as follow: 
Hex = LiJ (r— Rj). SiS (eq.43) 


Where S; and s represent the spin of the localized electrons and the delocalized carriers at Ri and r 


belonging to the magnetic ions and the semiconductor matrix, respectively. 


With the mean field approximation [95], the Hamiltonian exchange can be expressed in two forms: 
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For the symmetry s electrons in the conduction band and those localized on the magnetic ion: 
He, = —xNopa <S>s (eq.44) 
For the symmetry p electrons in the valence band and those localized on the magnetic ion: 
Hey = —XNoP < S>s (eq.45) 


Noaand Nofs are the exchange constants corresponding to the s-d and p-d exchange interactions, 
respectively. In the case of a positive exchange constant, the interaction between the spins is 


ferromagnetic else, it is antiferromagnetic. 


During the past ten years, these DMSs have been extensively studied, but there have been limited 
experimental investigations conducted, in contrast to the existing theoretical works available in the 
literature. Some notable examples include Cao.7sCro2sS and Cao.75Vo.25Se [96], Zni-xCrxS and Cdi-x.CrxS 
[97], Mgi-xVxSe [98], Beo.7sMno25X (X = S, Se, Te) [99], and Sro.7sTMo.2sS (TM is 3d transition metals) 
[100]. More remarkably, Fe-based DMSs have garnered considerable attention as ferromagnetic 
semiconductors because of their high Curie temperatures, low power consumption, and suitability for 
high-speed spin devices [101]. Researchers have been inspired by these intriguing characteristics and 
have dedicated efforts to unraveling the ferromagnetic mechanism in Fe-based DMSs. Several notable 
studies have contributed to our understanding of these materials. We mention, Fe-doped ZnO by Zhang 
et al. [102], Fe-doped CdS by Bourouis et al. [103], Fe-doped CdSe by Singh et al. [104], Fe-doped 
ZnTe by Mahmood et al. [105], Fe-doped ZnSe by Li et al. [106], and Fe-based CaS by Amari et al. 
[107]. The results of these studies showed that the materials exhibited their highest strength when 
incorporated with the Fe impurity, exhibited ferromagnetic behavior at room temperature, and allowed 


tuning of mechanical, magnetic, and optical properties through different Fe dopant concentrations. 


The co-doping process has shown even greater improvement in such DMS, as it has been 
demonstrated that the addition of other impurities besides the Fe atoms not only improves the physical 
properties but also introduces new properties different from those obtained by single doping can be 
achieved. An example of this is the addition of Fe alongside Mn in CdTe [108], which not only improves 
the magneto-optical properties but also changes the DMS type from n-type to p-type. Another example 
is the co-doping of CdS with (Cu, Fe) [109], which results in a transformation of its character from 
metallic to semi-metallic in nature with an increase in magnetic moment value, making it useful for 
spintronic applications. Furthermore, unlike other transition metals such as Cr and Co, the addition of 
Fe to CdxHg: - xSe [110] results in a significant change in thermoelectric properties, which could be 


manipulated by adjusting the composition. 
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On the other hand, the presence of TM clustering, the complexity of controlling stoichiometry, and 
the challenges associated with d-d transitions make the use of both TM- co-dopants difficult [111]. 
Nowadays, co-doping with other substituents such as non-magnetic elements, called dO, has emerged as 
a viable approach to improve and tailor the properties of DMS materials and offers exciting opportunities 
for advancements in various technological applications. Further discussion of this type of material can 


be found in the following section. 


1.4.2.2 d° vs. Conventional Magnetism 


d° magnets include a category of materials that lack magnetic ions with d or f orbitals and should 
theoretically not exhibit ferromagnetic properties. However, these materials often exhibit 
ferromagnetism with a Tc of over 300 K [112]. In this fascinating category of materials, the first 
materials discovered were doped non-magnetic oxides, thin films of HfO2, non-stoichiometric CaBe and 
irradiated graphite, with observation revealing that these materials exhibit small ferromagnetic moments 
and possess Curie points well above room temperature despite the absence magnetic atoms [113]. The 
magnetism in these materials has been attributed to the interaction between localized magnetic moments 
and traveling charge carriers. Even though d’ elements do not have partially filled d or f electron shells, 


they can still exhibit ferromagnetism due to the double exchange mechanism discussed previously. 


In some cases, the formation of resonant donor levels in d° ferromagnets can also contribute to their 
magnetic properties. These resonant donor values can be controlled by varying the composition of the 
solid solution, allowing tuning of the magnetic behavior [114]. From now on, a significant part of the 
work was devoted to the study of ferromagnetism, which arises from the incorporation of non-magnetic 
elements into wide bandgap semiconductors [115-117]. This increase in interest can be attributed to the 
remarkable effects of metal-free magnetism, which offers the potential for the further development of 
lightweight, biocompatible and environmentally friendly magnetic materials [118]. A widely discussed 
scenario for d’-based DMSs involves magnetic polarization of valence states by doping with sp-type 
impurities. To maximize HM ferromagnetism and improve structural and optoelectronic properties, for 
example, to increase luminance, facilitate charge injection, and achieve better energy band alignment, 
the most promising approach is to selectively dope the cation sites with the s-type, also known as p? 


elements from the first row of the periodic table, such as Lit, Na* and K* [119]. 


The alkali metals Li*, Na* and K* have one valence electron and are therefore electron donors when 
incorporated into a host material. This means that they provide the host material with additional charge 
carriers (electrons), effectively increasing its electrical conductivity. Doping with alkali metals can 
improve the electrical and thermal properties of the material, making it more suitable for various 
electronic and energy applications. In addition to changing the electrical conductivity, doping with alkali 


metals can also influence the optical properties of the material. For example, there may be changes in 
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the bandgap, allowing the material to absorb and emit light of different wavelengths. This property is 
crucial for the development of optoelectronic devices such as LEDs and solar cells. The monovalent 
elements Li*, Na* and K* are also in close proximity to the position of Sr in the periodic table. Due to 
their similar atomic size, they are well suited for incorporation into the rock-salt SrS lattice, which 
represents a high probability for co-doping with transition metals (TM) and has a significant influence 
on its properties. Surprisingly, investigations into alkali doping at the Sr site are relatively new. In an 
experimental study by Chang et al. Na*-doped SrS phosphor has been shown to produce high-quality 
solid-state light-emitting diodes (SS-LEDs) with improved emission intensity over the wavelength range 
of 430-700 nm [120]. Furthermore, the ferromagnetism in the compound (Sro.75 Ko.25 )24nGe2S60 was 
increased by replacing Sr** with K*, as theoretically predicted [121], resulting in a Tc of 300 K. Studies 
by Yang et al. [122] proved the effectiveness of cation doping in Li*-ion and Na*-ion batteries and 


showed promising results for large-scale electrical energy storage. 


While extensive experimental and theoretical studies have been carried out on the co-doping of 
alkali elements with transition metals in other matrices [123-126], the study of the SrS matrix provided 
with both co-dopings is completely missing. Therefore, it is crucial to study in depth the co-doping 


process within these elements and its effects on various SrS properties. 


After explaining the basic methods used in this dissertation to tune SrS bulk properties, we now 


provide a brief summary of dimensionality reduction, focusing on the 2D monolayer configuration. 


1.5 Low-Dimensional Materials 


Changing material properties on the way to the nanoscale is a field of intensive research. With its 
focus on materials and structures ranging in size from 1 to 100 nm, nanotechnology opens up a world 
of possibilities [127]. While the term “nanomaterial” may only conjure up thoughts of size, it 
encompasses much more. Nanotechnology deals with the manipulation of structures at the atomic and 
molecular levels and the development of materials that are precisely tailored to specific applications 
[128]. The advanced application of nanomaterials highlights that dimensions and shapes are critical 
aspects that have a profound impact on performance [129]. Nanomaterials are categorized based on their 
dimensions: zero-dimensional (OD), one-dimensional (1D), two-dimensional (2D), and _three- 
dimensional (3D) nanostructured materials (see Figure 14) [130]. This emphasis on precise design and 
careful structuring is motivated by the goal of overcoming challenges in various areas. Whether in 
wastewater treatment, energy production, conversion or storage, the targeted manipulation of these 


materials holds the potential to bring about transformative advances. 
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Figure 14. Taxonomy of nanoscale dimensions. (Source: Tallinn University of Technology). 


Among the diverse families of nanomaterials, 2D nanomaterials have emerged as a particularly 
fascinating class characterized by a variety of fascinating properties. In short, this category is discussed 


in the following section. 
1.5.1 Two-Dimensional (2D) Materials 


2D materials include a specific class of substances in which one dimension operates at the nanoscale, 
while the remaining two dimensions exceed this threshold. Over the past 16 years, intensive study of 
these materials has increased significantly, particularly those composed of layered compounds 
characterized by robust intralayer bonds and sensitive van der Waals forces between layers. As the 
electronics field experiences rapid growth mainly due to advances in silicon technology [131], scientists 
have been working hard to develop new materials that could revolutionize the technology of electronic 


devices while solving the complex problems of heat dissipation. 


In 2004, A.K. Geim and K.S. Novoselov of the University of Manchester reached a significant 
milestone with the successful isolation of graphene, the groundbreaking 2D material [132]. Their 
groundbreaking efforts earned them the Nobel Prize in 2010. Graphene, a crystalline carbon allotrope, 
has a zero bandgap and a two-dimensional atomic structure with a hexagonal arrangement at the atomic 
scale. In this structure, each carbon atom forms four bonds: three o bonds (sp” hybridized) with its 
neighboring atoms and one z bond oriented perpendicular to the plane. This fundamental structural unit 
serves as the basis for other carbon allotropes, including graphite, fullerene, nanotubes, and nanocones, 


making graphene the precursor to all carbon-based nanomaterials. 
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The discovery of graphene has sparked an explosion of research in the field and opened a new 
chapter in the study of other 2D materials. Scientists have created monolayers, bilayers, and materials 
with some layers that resemble graphene [133]. Along with other inventions, they have also created 
nanoribbons [134], single-walled and multi-walled nanotubes [135] and more. The unique ability of 
these materials to adjust their energy bandgap has enabled their use in a wide range of fields including 
electronics, optoelectronics, semiconductors, batteries, composite industries, solar energy, 


communications and more. 


The 2D materials space experienced even greater expansion with the introduction of Ti3C2T;, a 
notable member of the MXenes family. This groundbreaking development took place at Drexel 
University in 2011. Subsequently, similar to graphene, silicene formed as the first member of the 
MXenes group in 2012 [136]. This was the beginning of a series of discoveries that revealed various 
MXenes materials, which now include a diverse range of over 30 different 2D types. A quick 
examination of the dimensional database confirms the significant increase in research into 2D materials 


in recent years (Figure 15). 
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Figure 15. Annual publications on 2D nanomaterials over the past decade. 


The important surface-to-volume ratio of 2D materials as well as the exceptional occurrence of 
electron confinement enable precise modification of properties. High thermal conductivity, remarkable 
charge carrier mobility even at ambient temperature, and atomically thin properties are all found in 2D 
materials, which also enable simplified processing techniques. Due to their exceptional mobility, 2D 
materials are more sensitive to conductivity. This is because surface-generated supports respond 


dynamically to photoeffects. This intriguing property positions them at the forefront of high-gain 
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photodetector applications, spanning areas such as optical communications, optoelectronic devices, and 


biomedical imaging. 


Furthermore, their wide electromagnetic absorption range, ranging from infrared to ultraviolet 
[137] makes 2D materials central components in high-performance photonics and optics. Their ability 
to maintain low absorption (<10%) yet exhibit high conductivity makes them invaluable for transparent 
electronic applications, including solar cells and liquid crystal devices [138]. Beyond their remarkable 
electrical and photonic properties, 2D materials exhibit outstanding mechanical properties. These 
materials are characterized by their remarkable flexibility and have a breaking strength 200 times higher 


than that of steel [139], making them crucial for reinforcing polymers in composites. 


1.5.2 Synthesis Approaches for 2D Materials 


Two approaches are used in the synthesis of 2D materials: the top-down and the bottom-up approach 
[140]. A comprehensive overview of various top-down and bottom-up approaches is shown in Figure 
16. In the coming discussion, we will briefly examine these approaches and highlight their specific 


features. 


1.5.2.1 Top-Down Approach 


In the top-down approach, the process begins with bulk materials, layers of which are carefully 
stripped or peeled away using various techniques to obtain atomically thin 2D layers. A notable example 
of this approach is peeling off layers of graphite using adhesive tape, creating the ultra-thin monoatomic 
layer called graphene, which is made up of graphitic carbon atoms. In addition to the methods shown in 
Figure 16, the transition from 3D bulk structures to 2D structures is achieved through a variety of other 
strategies. These include, in particular, liquid phase peeling [141], ion intercalation/peeling [142] and 
physicochemical peeling [143]. Among these, chemical exfoliation stands out as a widely preferred top- 
down approach to producing materials with small dimensions. Several other top-down techniques are 
used, including mechanical milling [144] and lithography, including techniques such as ion beam 
lithography, electron beam lithography and photolithography [145]. In parallel to these methods, 


approaches such as sputtering [146], the arc discharge process [147] and laser ablation [148] are used. 


1.5.2.2 Bottom-up Approach 


In the bottom-up approach, atoms are cleverly joined together, similar to sewing, to create 
atomically thin 2D material structures. Numerous bottom-up techniques are used, including molecular 
beam epitaxy, chemical vapor deposition (CVD), physical vapor deposition (PVD), hydrothermal and 


solvothermal methods, and template-based methods [149]. Bottom-up approaches are known for their 
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ability to produce high-quality 2D structures and find application across a spectrum of 2D materials. In 
particular, they play an important role in the design of two-dimensional metal-organic frameworks 


(MOFs), coordination polymers (CPs), and covalent organic frameworks (COFs) [150]. 
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Figure 16. Widely utilized top-down and bottom-up approaches to synthesize the 2D materials. After 
Ref [151]. 


1.5.3 Categorization of 2D Materials 


Due to their finely tunable energy bandgap and remarkable intrinsic properties, 2D materials have 
become key players in a spectrum of innovative fields, including electronics, optoelectronics, 
semiconductors, spintronics, spin caloritronics, batteries, composite industries, solar energy utilization 
and communication systems, and a range of other applications [152]. Numerous subclasses within this 
material category, such as silicene [153], germanene [154], hexagonal boron nitride (h-BN) [155], 
stanene [156], phosphorene [157] and single-layer transition metal dichalcogenides (SL-TMDs) [158] 


were included Focus of intensive and extensive research efforts. 


These 2D materials form the elementary building blocks for critical components in low-dimensional 
devices. Based on this, we can divide 2D materials into three main families accordingly: the graphene 
family, a large class that includes graphene, -BN, BCN, fluorographene, graphene oxide, and more; the 
2D chalcogenide family, an ensemble that includes MoS2, WSe2, ZrS2, NbSe2, Bi2Ses, etc.; and finally 
the family of 2D oxides, a broad category that includes hydroxides such as Eu(OH)2 and layered Cu 


oxides, MoO3, WO3 and others, making a fascinating representation as shown in Figure 17. 
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Figure 17. Schematic illustration of the different 2D nanomaterials families. After Ref [151]. 


The group of 2D chalcogenides containing group VIA elements is currently one of the most 
important groups and offers a wide range of materials, crystal structures and properties. This has 
attracted increasing attention due to their large bandgap, abundance, and special properties, including 
high mobility and efficient absorption in the visible and ultraviolet regions. 2D chalcogenides, either in 
their monolayers, bilayers, or multilayers, can be further categorized based on their chemistry and 
stoichiometry, as shown in Figure 18. These materials are made up of atomic planes, which are self- 
contained units with no dangling bonds present on surfaces. This special property facilitates the isolation 
of individual layers from the main material and shows the remarkable potential of these compounds 
[159]. 
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Figure 18. Categorization of 2D chalcogenides. After Ref [159]. 
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A significant subgroup within the chalcogenide family are the AX monochalcogenide monolayers, 
essentially the II-VI group, which have also attracted considerable attention from researchers in recent 
years. This lattice arrangement has a honeycomb structure consisting of A-X bonds (where A means a 
cationic element from the group of If elements and X is the chalcogen elements (S, Se and Te), and is 


endowed with the Dn point group symmetry. Each atom forms a vertex through sp hybridization. 


In Figure 19 we show a visualization of the hexagonal structure of SrS, a prototype of the II-VI 


family. This symmetry is similar to other known monolayers, including graphene, silicon and BN. 


Figure 19. a) The Unit cell of the monolayer SrS, which adopts a hexagonal structure. Green spheres 
represent the Sr atoms, while orange spheres depict the chalcogen S atoms. The image was generated 
using VESTA software [76]. b) The first brillouin zone of a hexagonal crystal, with the irreducible 


wedge identified by their common names, and the reciprocal base vectors labeled as bi, and bo. 


Similar to their traditional 3D II-VI compounds, density functional theory (DFT) was also used to 
predict the novel 2D counterparts. For example, the 2D monolayers BeO, MgO, CaO, ZnO, CdO, CaS, 
SrS, SrSe, BaTe and HgTe were found to have strong dynamic and thermal stability and are suitable for 
use in UV optoelectronics [160]. The CdO monolayer has also shown promising optoelectronic 
capabilities suitable for visible light photocatalysis [161]. Recently, Abdullah et al. [162] showed that 
the combination of Be with O results in different thermal and optical properties, including high heat 
capacity and an active optical response observed in the UV. Interestingly, SrS and SrSe showed good 
thermal and dynamic stability and exhibited shallow valence bands that became highly spin polarized 
when doped with holes, making them excellent options for realizing Stoner ferromagnetism [160]. More 
recently, Yari et al. [163] found that doping the SrS monolayer with Cr induces a pronounced HMF 
character and increases the thermoelectric quality factor (ZT) to 0.9 at temperatures above 500 K, 


making this monolayer suitable for spintronic and spin caloritronic nanodevices. 
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Reviewing the existing literature, the need to develop a diluted SrS monolayer-based magnetic 
semiconductor appears to be extremely promising. If we focus our attention on SrS as the original 
compound, the study of dimensionality reduction becomes extremely compelling. The aim of this study 
is to understand the unique influence of dimensionality reduction on various properties of SrS, 


particularly in the context of Fe doping compared to co-doping process. 


Therefore, a bibliographic review of bulk SrS, monolayer SrS and their doped or co-doped 


structures is required. 
1.6 An Overview of Bibliography 


In the last decade, spintronics and its derivatives, including spin caloritronics, have come to the 
forefront of research and include both bulk (3D) and two-dimensional (2D) structures. In particular, the 
group of II-VI semiconductors based on transition metals (TM) proves to be a central component. With 
the ability to cover the properties of metalloids and semiconductors, these materials have a wide range 
of properties and applications by simply changing their composition. Their characterization, structural 
phase transition, synthesis techniques, magnetic behavior, electrical and optical responses, 
thermoelectric characteristics, and a variety of other properties have all been the subject of extensive 


research [164, 165]. 


In the following discourse, we discuss some of the groundbreaking work carried out on the TM- 


based II-VI semiconductor SrS, highlighting both the 3D and 2D forms. 


e The binary compound SrS has been subject to thorough investigation, both experimentally and 
theoretically, in which it is confirmed that this compound shows p-type semiconductor 
properties with an indirect bandgap in both its 3D bulk structure and its 2D monolayer 


counterpart [73, 75, 82, 166-172]. 


e Xiao et al. [173] have examined the thermoelectric characteristics of bulk SrS, while Yari et al. 
[174] and Rajput et al. [175] have directed their attention towards examinining the 
thermoelectric quality at the monolayer level. The collective findings underscored the promising 


potential of SrS for application in thermoelectric devices, with the ZT value approaching unity. 


e The electronic structure, half-metallicity, and ferromagnetic properties of doped-SrS have been 
subject to theoretical investigations by Bourega et al. [176], Doumi et al. [64, 65], and Hamidane 
et al. [177] utilizing an ab-initio study. The combined results have demonstrated a large half- 
metallic bandgap and a stabilization of the ferromagnetism through exchange interactions 
between the host and magnetic dopants. These results underline how well these compounds 


might work as spintronic device materials. 
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e Hamidane et al. [66] conducted ab-initio calculations on the Srj..MnxS bulk systems 
(0.25<x<0.75). Their study revealed that the introduction of Mn into the host SrS induced 
alterations in the band structure and density of states. Interestingly, these changes enhanced the 
optical characteristics with increasing Mn concentrations, induced a significant magnetic 
moment, and preserved the material's semiconducting nature. 

e Inhis research, Hoat [100] explored the electro-magnetic properties of SrS doped with first-row 
transition metals, maintaining a constant concentration of x = 0.25. His findings indicated that 
among all materials, only five, including Sro7sFeo25S, could potentially serve as candidates for 


spintronics, primarily due to the presence of a considerable half-metallic bandgap. 


e The hole doping at the anion site (with P and As) and cation site (with Al and Ga) of SrS 
monolayer was first investigated by Lin et al. [172], where they observed that the doping 


strategy resulted in an increase in the magnetic moment and the spin-polarization. 


e Inarecent study, Yari et al. [164] demonstrated the existence of magnetism in SrS monolayer 
through hole doping with Cr. Through a comprehensive exploration of its electronic, magneto- 
optical, and thermoelectric properties via a DFT approach, he confirmed the suitability of this 
monolayer for applications in spintronics, spin caloritronics, and spin-based optoelectronics 


devices. 


Besides these, diverse nanostructures involving doping magnetic elements into SrS have been grown 
experimentally using different techniques such as sol-gel route [178], electrochemical deposition 


[179], solid-state diffusion [180], and others [181-183]. 


1.7 Simulation Studies and Their Scope 


Broadly speaking, the study of structures involves two main approaches: theoretical and 
experimental. However, the advancement of computing capabilities has led to a distinct path known as 
Computational Condensed Matter Physics. This innovative approach acts as a mediator and connects 
the areas of theory and experiment. The development of computational resources has greatly facilitated 
deeper exploration of matter across different dimensions and length scales. Computer simulations are 
becoming increasingly common in condensed matter physics, both at the macroscopic and mesoscopic 
scales. Macroscopic systems are systems that are both measurable and observable, while mesoscopic 


systems deal with systems ranging from nanometers to micrometers. 
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Computer simulations are sometimes referred to as computer experiments because they share many 
similarities with laboratory experiments. The beginning of a computer simulation involves constructing 
an idealized model of a particular physical system. An algorithm or method is then defined to implement 
this model on a computing platform. By executing a computer program, the physical system is simulated 
and the basic ideas of the computer experiment are summarized. This digital experiment serves as a link 


between theoretical models and laboratory studies. 


1.7.1. Possible Simulation Approaches 


A variety of simulation methods have proven invaluable in gaining insight into the field of 


condensed matter physics, generally categorized as follows: 


e Classical Molecular Dynamics: This method employs classical (Newtonian) mechanics [184], 
enabling the dynamic representation of systems. 

e Quantum Monte Carlo: Characterized by a stochastic approach, Quantum Monte Carlo stands 
as an almost exact method for determining the many-body wavefunction [185], providing a 
robust avenue for exploration. 

e Ab-initio Method: Grounded in computer simulation, this approach adresses the quantum 


mechanical Schrédinger equation [186], providing a complete solution. 


Based on independent electronic structure computations, classical molecular dynamics is a 
powerful and well-established technique for studying many-body condensed matter systems [184]. 
The interatomic potential approximation is a fundamental question in any molecular dynamics 


scheme. Typically, the molecular dynamics method entails predetermining these potentials. 


Quantum Monte Carlo (QMC) is revealed as a stochastic technique for calculating the ground 
state parameters of quantum systems by solving the Schrédinger wave equation. The QMC 
technique has shown promise in a number of fields with Schrédinger-like Hamiltonians, including 
solids, quantum liquids, nuclear matter, spin systems, and ab-initio quantum chemistry. Using 
significant sampling techniques is a key component of QMC approaches. This involves using a trial 
vector, based on "zero-variance property", which indicates that faster convergence toward the exact 
eigenvector is caused by a tighter alignment between the trial vector and the exact one, resulting in 


less statistical fluctuations [185]. 


Finally, ab-initio methods utilize first-principles to compute material properties by numerically 
solving the quantum mechanical Schrédinger equation [186]. These methods serve as primary tools 
for conducting research in condensed matter physics, materials science, quantum chemistry, and 


molecular chemistry. The level of accuracy predominantly depends on how accurately one can solve 
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the electronic problem. Among the most accurate based-methods, the Density Functional Theory 
(DFT). Density Functional Theory (DFT)-based methods include techniques that directly solve the 


eigenvalues and eigenfunctions of the molecular electronic Hamiltonian (H) using basis set methods. 


In Chapter 2, a detailed study of the DFT approach and the techniques used in this dissertation 
is presented. The application of DFT has proven crucial for studying a wide range of ground state 
properties, spanning both bulk and low-dimensional structures. Therefore, in this context, we 


adopted a theoretical approach to thoroughly investigate the systems selected for our study. 


1.8 Conclusion 


Finally, this chapter provides an introduction to II-VI bulk 3D materials and their corresponding 2D 
monolayer derivatives, followed by a detailed literature review of their properties and applications. Due 
to their remarkable electronic properties, characteristic structures, dimensionality, high carrier mobility, 
and tunable band gaps, these materials have attracted great interest as potential alternatives in various 
fields. Advances in storage technology and the constant pursuit of environmentally friendly solution 
materials provide the scientific community with numerous new opportunities to explore phenomena, 
concepts and technologies in the fields of spintronics, spin caloritronics and spin-based optoelectronics. 
Dilute magnetic semiconductors (DMS), obtained by incorporating transition metals into I-VI 
semiconductors, hold enormous potential for development into small-sized, low-power, and high- 
capacity memory systems. Through its comprehensive discussions, this chapter contributed significantly 
to the selection of suitable host semiconductors and dopants for our dissertation. We specifically chose 
to focus on Fe-doped SrS, using two different strategies to improve its properties: co-doping in the 3D 
form and dimensionality reduction to the 2D monolayer form. These approaches aim to improve a wide 
range of properties including structural, electronic, magnetic, mechanical, optical and thermoelectric 


properties. 
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Chapter 2 


Theory and Computational Techniques 


his chapter presents a fundamental theoretical formulation to address the challenges of the 
many-body problem through the application of density functional theory (DFT) and a 
computational formulation for ab-initio calculations. It also describes the software and 


methods used in the calculations carried out as part of our study for this dissertation. 


2.1 Introduction 


The host material investigated in this work is a semiconductor. The study of a semiconductor or 
any solid begins with the study of an ideal crystal at the atomic level at a temperature of 0°C. This crystal 
is made up of a series of atoms (or ions) that occupy precise positions and repeat periodically to form 
the material. The study is based on the interaction between electrons themselves and with nuclei. Due 
to the nature of particles (electrons and nuclei), their study requires the use of quantum mechanics, which 
requires solving the Schrédinger equation to determine the total energy of the system. All properties of 
materials can be determined by using appropriate computational tools to solve their quantum mechanical 
problem. This theory governs the electronic structure, which is responsible for other solid-state 
properties such as optical, electrical, magnetic, mechanical properties, etc. Electrons and nuclei form a 
highly complex N-body system, which makes solving the Schrédinger equation extremely difficult or 
even impossible. Various methods have been proposed to solve this problem. The method that has 
proven to be extremely successful and is the most widely used is density functional theory (DFT). With 
the introduction of density functional theory (DFT), it becomes possible to describe the particle system 
(nuclei and electrons) on the basis of fundamental data: lattice parameters and atomic numbers of the 


elements. 


This theory allows us to gain insights into the composition of materials and predict their potential 
applications, thereby helping experimentalists develop novel devices. The study of the fundamental 


properties of atoms, molecules, solids and systems with reduced dimensions in their ground states forms 
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a central focus in the field of condensed matter physics. The fundamental understanding of the electronic 
structure of materials depends on solving the Schrédinger wave equation. This equation effectively 
illustrates the behavior of uncomplicated systems (such as the hydrogen atom problem), and numerical 
solutions are only possible for a limited number of atoms and molecules. However, when complex 
systems are involved, solving this equation and then estimating physical properties becomes extremely 
challenging. Although we know how to solve the problem, we still lack powerful tools to find the 
solution. Over the past century, there have been significant changes in the way we solve the Schrédinger 
wave equation for systems with many particles. The first breakthrough came in 1928 when Hartree 
approximated the wave function for many particles and gave an exact energy value for the hydrogen 
atom (-13.6 eV) [1]. Thereafter, efforts were made to find a better wave function using Slater 
determinants, but dealing with large systems remained difficult due to the need for extensive 
computational resources. The real progress came when three-dimensional electron density theories and 
energy function methods were introduced. This approach not only reduced computational effort but also 


proved remarkably accurate. 


This chapter examines the development of this theory and provides a detailed explanation of the 


Kohn-Sham approach and the approximations used for exchange and correlation. 


2.2 Theoretical Formulation 


2.2.1 Many Body Problem within Schrédinger Wave Equation 


Here we introduce the basics of many-body quantum mechanics. We begin by explaining the basic 
principles of quantum mechanics and the properties of wave functions. These ideas will later help us 
reformulate the problem in terms of electron density. We start with the Born-Oppenheimer 
approximation and other simplifications to solve the Schrédinger wave equation for systems with many 
particles. Finally, we examine how density functional theory (DFT) is formulated and applied to study 
various properties of materials in different dimensions, from bulk materials to low-dimensional 
structures. In their solid state, materials are composed of electrons, which have tiny mass and negative 
charge, and nuclei, which carry significant mass and positive charge compared to electrons. The 
macroscopic properties of materials are determined solely by the positions of these electrons and nuclei. 
As a result, electrons in materials exhibit quantum behavior and their interactions are closely linked by 
quantum mechanics. The complexity of the quantum mechanical system under consideration increases 


significantly as the number of electrons increases, resulting in the so-called many-body problem. 


To put it more precisely: If there are N nuclei/electrons, there are a total of N+ZN interacting 
particles. These particles come into contact with each other through electromagnetic forces. Due to the 


tiny size of these interacting components, the system under study must be viewed from a quantum 
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mechanical perspective. Traditionally, the properties of a quantum mechanical system are determined 


by solving the time-independent Schrédinger wave equation, which is given as follows [2]: 


Aw(r,, To, ... Ry, Ro...) = EW(r,, 72, ... Ry, Ro, --) (eq.1) 


Here, the variables 77, R, represent the position vectors of electrons and nuclei, respectively, H 
denotes the many-body Hamiltonian, EF is the eigenvalue of energy and w represents the eigen 


wave function. 


The Hamiltonian (H) in connection with many-body problems is based exclusively on electrostatic 
interactions (electron-electron, electron-nucleus, nucleus-nucleus) and kinetic energies. This essential 
quantity can be described with the following equation: 


Z Ze" 
IR;-Ry| 


h? D2 h? 
_ 2m, de Vi = ie 2m, = Sy (eq.2) 


Ir,-Ry| i + 3D = + ys 


The indices i andj are used for electrons and J and J for nuclei. On the left side of the equation provided, 
the first two terms represent the kinetic energy of electrons and nuclei, respectively. The remaining three 
terms characterize the energies associated with interactions: electron-nucleus, electron-electron, and 
internuclear Coulomb interactions. If we assume atomic units denoted as (h? = e? = m = 4mt€9 = 


1), the Hamiltonian is further simplified to: 


ZZ; 


oe 1.S> 1 
H =-> iV? - ye v2 - dil ar + > Lieje= = = Hae De = (eq.3) 
In terms of the operators, Eq.3 can be formulated as follows: 
H=T. +Ty + Ven + Veo +Vyw (eq.4) 


Our goal is to investigate the many-body Schrédinger wave equation for a complex system using various 
techniques and approximations. To simplify this Hamiltonian, we use the Born-Oppenheimer 


approximation, which separates electronic and nuclear motions. 
2.2.1.1 The Born-Oppenheimer Approximation (1927) 


The forces on electrons and nuclei due to their electrical charge are comparably strong, leading to 
similar changes in their momenta due to these forces. This means that their impulses are comparable in 
size. Due to the significant difference in mass between nuclei and electrons, the velocities of the nuclei 
are significantly slower than those of the electrons. Consequently, we can assume that the nuclei are 


essentially stationary. This allows us to first solve the Hamiltonian for the electronic ground state and 
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determine the energy of the system in this configuration before dealing with nuclear motion. This 
separation between electronic and nuclear motion is called the Born-Oppenheimer (BO) approximation 
[3]. The separation between electronic and nuclear motion allows us to represent the complete wave 


function as a multiplication of electronic and nuclear wave functions, like this: 
w@,R) = O@R)x(R) (eq-5) 


Where 7(R) represents a nuclear wave function and 9(F, R) represents an electronic wave function 


that corresponds to certain nuclear positions. As a result, the Schrodinger wave equation for an electronic 


Hamiltonian takes the following form: 
A.(R)9@,R) = (Te + Vee F) + Ven @R) PCR) = En(R)OG,R) (eq.6) 


This produces a series of normalized eigenfunctions (7, R) and eigenvalues Ey (R) that depend 


on the kernel positions R. For each solution, there is a kernel eigenvalue equation: 


(Ty + Dyn @) + €n(R)) X(R)=E X(R) (eq.7) 


Since many of the desired properties can be derived from the ground state wave function and the energy 
of the system, our focus is on finding a technique to determine the ground state corresponding to the 
lowest energy state of the system. The variational principle provides a method of reaching the ground 


state by minimizing the total energy of the system. 
2.2.1.2 Variational Principle: Ground State of the System 


Using the variational principle, the energy calculation derived from an approximate wave function 
y(T) acts as an upper bound on the true ground state energy Eo. Let us consider a system that exists in 


a state |y >. The energy is estimated by evaluating the expectation value of the Hamiltonian and is 


represented as follows: 


EQ) = by (eq.8) 


By using the estimated value of the energy expectation (E) in the state |y >, it becomes possible to 


approximate the eigenvalue and eigenfunction of energy using the variational principle. 


SE(w) = 0 (eq.9) 
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The variational principle helps approximate the ground state energy of a system. Therefore it is 
important for ground state calculations. However, the techniques for implementing this principle can 
vary depending on the approach used to estimate the wave function and associated energy for the system. 
Numerous approximations have been developed to find the optimal wave function that corresponds to 
the antisymmetry property and includes various electron interactions within the system. To estimate the 
desired wave function, the Hartree and Hartree-Fock approximations are used, which will be discussed 


in more detail below. 


2.2.1.3 Independent Electrons Approximations 
a) Hartree Approach (1928) 


The Born-Oppenheimer approximation simplified the complexity of the Schrédinger wave 
equation, but the challenge remained due to electron-electron interactions. Hartree addressed this 
complexity by treating electrons as separate entities, resulting in the Hartree approximation [1]. In this 
approach, the overall wave function is represented as a product of individual one-electron wave 


functions. Consequently, 


o(F, R) = Q—i (7, R)9, (72, R) oo Qu (Ty, R) (eq.10) 
M =. 
= || Px (Fx, R) (eq.11) 
k=1 


This is known as the Hartree product. Here the functions @, (Ty, R) are also called orbitals and 
must satisfy the condition of being orthonormal to each other, which can be expressed as 
f dr 9.1) @).(7) = Sx. This approximation does not satisfy the antisymmetry principle, also 
known as Pauli's exclusion principle [4]. This principle states that the wave function describing fermions 
must exhibit antisymmetry when exchanging any set of spatial and spin coordinates. Therefore, Hartree 


theory requires adjustments, a role that Hartree-Fock theory takes on. 
b) Hartree-Fock Approach (1930) 


The emergence of the Pauli exclusion principle led to the development of the Hartree-Fock method, 
which allowed the representation of the overall Hamiltonian for an N-electron system (2) as the sum of 
individual electron Hamiltonians (Hi), meaning = )), H; . Furthermore, the overall wave function is 
expressed as a Slater determinant, which consists of individual electron wave functions. In the case of 


an N-electron system, the wave function can be formulated as follows: 
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1%) X2(%2) XN) 
yp = a 1%) X2(%2) = tn Gn) (eq.12) 
X1%1)  X2(%2) XN) 
The prefactor of a guarantees the fulfillment of the normalization requirement for the wave 


function @. The determinant represented in Eq. 12 is recognized as a Slater determinant [5]. In the 
equations provided previously, we switched our notation from using spatial orbitals (@ (r)) to 
introducing spin orbitals (7(x)). Therefore, a spin orbital is essentially the result of multiplying a spatial 
orbital by a spin function, represented as (x) = @ (7) o. Here x = {r, o} denotes the combined set of 
space-spin coordinates, where r stands for the space coordinate and o stands for a spin coordinate. This 
spin coordinate can take one of two possible values: 61() or 62(|). This approximation gives the full 


energy as the Hartree-Fock energy (EHF), which is represented as follows: 


This approximation gives the full energy as the Hartree-Fock energy (EHF), which is represented 


as follows: 


1 1 


1 1 
Enr = (xilh|xi) + zi (xix; xix) _ zi (nix; xix) (eq.13) 


‘ij ‘ij 


In the context of the first term, # represents the one-electron operator and can be formulated as 


follows: 
F 1 
h(i) = Te +Vee = — 5Vi = yap (eq.14) 
The second term refers to a two-electron system and its integral form is as follows: 


1 


XiXxj Fj Xixj 


7 yy - ax'Xi (x) Xj Oxi (x) xj(x') 
lr—7'| 
ij 


-(x)I2 |v; (xr) |2 
= 3, f dxda! MOPED 


= 3) 


[FF 
= Yili (eq.15) 


Where Jj; is called the Coulomb integral. The third term also refers to a two-electron system and its 


integral form is as follows: 


1 aN aya 4169) 4108) 413) 1c @) 
XiXxj Ty xjxi) =) xdx aaa = aan 
uj 
= diy Ky (eq.16) 
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Where Kj; is identified as Exchange Integral. 


This approximation misses the Coulomb correlation, causing the total electronic energy to deviate 
from the exact solution. Consequently, the Hartree-Fock energy constantly exceeds the exact energy, 
and this deviation is called the correlation energy [6]. In addition to the Coulomb correlation, solving 
the Schrédinger wave equation presents another major challenge, known as the 3N variable problem. 
These problems of correlation energy and complexity due to 3N variables have been addressed by an 


innovative method called density functional theory (DFT), which is examined in the next subsection. 
2.2.1.4 Early Density Functional Theories: The Thomas-Fermi Model 


First, in 1927, Thomas and Fermi (TF) suggested that atoms could be imagined as uniformly 
distributed clouds of negatively charged electrons surrounding nuclei in a six-dimensional phase space 
that includes momentum and coordinates. This represented a significant simplification of the 
complicated many-body problem. It is instructive to explore the basic concepts of the TF approximation 
before delving into a more precise theory, namely density functional theory (DFT). The Thomas Fermi 
method [7] takes a different approach based on the electron density within the system. This method 
assumes that electrons move in an external potential. The main goal is to calculate the potential V(7")and 
the electron density (7). Electrons are assumed to be uncorrelated and their kinetic energy is estimated 
using a local approximation based on the behavior of free electrons. The potential can be found by 
solving the Poisson equation, and the requirement for a constant chemical potential leads to the Thomas- 
Fermi equation for n(7). Although this method is not successful for real systems, it serves as a prototype 


for density functional theory because it focuses on the electron density n(7), which is expressed as: 
n(r) = Nf dr, eeacane! ax dry|wcr, T 2, wennee aan ple (eq.17) 


The resulting energy functional includes classical expressions for nucleus-nucleus and electron-electron 


potentials, which were represented as follows: 
2 5 = 
Erpln@)] = 3 (3n?)s [ ns@ar — z [Par +5 If nO td, (eq.18) 


This model provides a good solution to the Schrédinger equation, but a clear relationship between 
n(r) and w(%) is required. A refined version of the Thomas-Fermi model was introduced by 
incorporating exchange effects via an exchange functional. However, this theory still gave inaccurate 
results for systems with many particles. The problem was finally solved when Hohenberg and Kohn 
showed that knowledge of the ground state electron density n(7) for any electronic system, independent 


of interactions, uniquely defines the system. 
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2.2.1.5 Density Functional Theory: An Ab-initio Approach (1964) 


This model provides a good solution to the Schrédinger equation, but a clear relationship between 
n(r) and w(r) is required. A refined version of the Thomas-Fermi model was introduced by 
incorporating exchange effects via an exchange functional. However, this theory still gave inaccurate 
results for systems with many particles. The problem was finally solved when Hohenberg and Kohn 
showed that knowledge of the ground state electron density n(7) for any electronic system, independent 


of interactions, uniquely defines the system. 


e Theorem I: In a system of interacting particles under an external potential Ve,,(7), the ground 
state density N9(7) uniquely determines the external potential. 

e Theorem II: Using the electron density n(7), a universal functional for the energy E[n] can be 
set up that is applicable to any external potential Vex. For a given Vex, the ground state energy 
of the system is the minimum value of the energy functional, and the corresponding electron 
density n(7), which minimizes this functional, is the exact ground state density. Thus, the total 


energy of the system can be expressed as: 
E[n] = F[n] + f 2? Vee @n@) (eq.19) 


The function F[n] = T[n]+Ve.e[n] includes both the kinetic energy and all electron-electron 
interactions. This functional is universal and independent of external potential. Therefore, ideally it 
should remain consistent across all systems. However, the theorems do not provide a method to 
determine the exact structure of this functional. Therefore, practical applications require approximations 


where the functional must be approximated for feasible calculations. 
2.2.1.6 The Kohn-Sham Formalism (1965) 


The Kohn-Sham equation [10] represents the Schrédinger wave equation for a fictitious system 
consisting of non-interacting particles, typically electrons. This constructed system yields the same 
density as any given set of interacting particles. The Kohn-Sham equation is characterized by a locally 
effective (imaginary) external potential in which the non-interacting particles move. This potential is 
usually referred to as Vxs(1") or Veg(7") and is called the Kohn-Sham potential. Since the particles in the 
Kohn-Sham system are non-interacting fermions, the Kohn-Sham wave function becomes a single Slater 
determinant formed from a collection of orbitals. The practical implementation of DFT was proposed 
by Kohn and Sham. The concept of the Kohn-Sham approach is to replace the interacting many-body 
system with a complementary system of non-interacting particles that have the same ground state. 


Consequently, the total energy functional can be expressed as: 
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E[n] = f d°F Vexe(F) NG) +7 ,(n) +3 f Fa? “AI se, (eq.20) 
Where Vex; is the external potential; 7; is the kinetic energy term of the hypothetical non-interacting 
electrons; The third term is the classical electrostatic energy (Hartree) of the electrons and all the many 
body effects are summarized in the exchange correlation energy term Exe. By minimizing the Kohn- 
Sham energy function with respect to the electron density n(7), a Kohn-Sham equation is derived that 


is similar to the Schrédinger equation: 
— — 1 — — = 
Hs@)piF) = |-5 V2 + Vis@ | vi® = ei (eq.21) 


The potential V5, in the above equation is an effective potential that includes the external potential 


Vex, the Hartree potential Vy, which is represented by the integral of f d?7’ ar and the exchange 
5E 
correlation potential Vy¢ = Soe where electrons move independently of each other. w,; represents 


the eigenfunctions corresponding to the eigenvalues €;. Consequently, the Kohn-Sham Hamiltonian in 


Eq.21 can be expressed as: 
= 1 — — — 
Hxs(F) = —5V? + Vext(T) + Vn) + Vc) (eq.22) 


This Hxg is related to the functional derivative of energy as follows: 


OE 


we HxsWi(r) (eq.23) 


The Kohn-Sham equations are theoretically precise. However, the functional form of Exc[n(7)] 


remains unknown and requires further approximations. The total energy contribution is made up of four 
terms (Eq. 22): the kinetic energy (-5V?), and the called the effective potential (Vore(1)), that 
englobes the classical Coulomb potential (Vex;(7)), the Hartree potential (Vq(1)) and the exchange 


correlation potential (Vy¢(7)). The exchange correlation potential (Vy¢(7)) includes a variety of many- 


body and quantum effects. Let us briefly discuss these terms one by one. 
e Kinetic Energy 


Kohn and Sham introduced a series of orbitals from which electron density can be calculated. It is 
important to note that these Kohn-Sham orbitals generally do not directly represent actual electron 
density. However, the connection between the Kohn-Sham orbitals and real electronic wave functions 
is that both give the same charge density. To calculate the kinetic energy term, the Kohn-Sham orbitals 


are used, as shown below: 
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T,[n] = DM f dF. Wi@(-[V WM (eq.24) 


e Coulomb Potential: Substituted with a Pseudopotential 


The V,,x¢(%) term characterizes the interaction of the electron with various atomic nuclei present, 
which essentially boils down to the Coulomb potential. To deal with this aspect, a technique is used in 
which the Coulomb potential is replaced by a pseudopotential [10]. It is generally accepted that core 
electrons are minimally involved in bond formation and most material properties depend primarily on 
valence electrons. Therefore, it is more practical to focus on valence electrons. These valence electrons 
experience the Coulomb potential, which is exerted by an inert ionic nucleus (pseudonucleus, consisting 
of core and core electrons). Consequently, the system is simplified to handle fewer electrons, as shown 
in Figure 1(a). Essentially, a pseudopotential serves as an approximation for the total nuclear potential. 
The effectiveness of a pseudopotential depends on how closely it reproduces the actual results. When 
creating a pseudopotential, it is necessary to ensure that the all-electron wave function agrees with the 


pseudowave function within a certain limiting radius, as shown in Figure 1(b). 


Nucleus 
Core electrons S / 
Valence electrons — 


(a) Pseudopotential model: The (b) Schematic representation of an atomic 
movement of external electrons occurs all-electron wave function (solid line) 
within a_ stable configuration of next to the corresponding atomic 
chemically inert ionic nuclei that pseudo-wave function (dashed line), 
include both core and core electrons accompanied by its corresponding 
[11]. external Coulomb potential and 


pseudopotential [11]. 


Figure 1. Pseudo potential approach — pseudo wave function and pseudo potential representation. 
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e Hartree Potential 


The Hartree potential describes the Coulomb repulsion experienced by the electron (at position 7), 
as shown in the Kohn-Sham equations. This repulsion is determined by the total charge density 
encompassing all electrons (positioned at 7’) within the given scenario. The expression for the Hartree 
potential is as follows: 

Va=f apse. (eq.25) 

To solve the KS equation we need the Hartree potential, which is closely linked to the charge 
density of the system. To determine the charge density, we rely on single-particle wave functions, which 
can only be obtained by solving the KS equations. This leads us to a challenge of solving a set of self- 
consistent equations. This process involves one iteration starting with an initial set of single-particle 
wave functions to derive the effective potential. The solution obtained using this method is called a self- 


consistent solution. The approach to solving this complicated equation is systematically shown in Figure 
2. 


Initial guess of density 


no(r) 


Generate effective potential 


Ver. (r) = Voxt.(T) t Vitar. (T) t Vec(Y) 


Solve Kohn-Sham equations 


Axsvbx(r) = exe (r) 


Calculate new density 


Ne 
n(r) = >> ba(r)? 
k 


Reached self-consistency? 


Done! 


Figure 2. Flowchart illustrating the various steps within the Kohn-Sham iterative loop in Density 
Functional Theory (DFT) codes. 
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e Exchange-Correlation Potential 


The most important potential term in KS equations is Vxc, which is used to account for exchange 
and correlation effects. In the following sections, we first deal with the elucidation of exchange and 
correlation. We then move on to discussing various approximations that are essential for estimating 


exchange and correlation. 
2.2.1.7 Exchange and Correlation 
2.2.1.7.1 Exchange-correlation Hole 


The principle of orbital antisymmetry dictates that electrons with the same spin must inhabit 
different orthogonal orbitals, requiring spatial separation between them. This results in a reduced 
electron density called an exchange hole, which results in reduced repulsion and consequently reduced 
system energy. This phenomenon causes the electrons to move away from each other, effectively 
eliminating the self-interaction inherent in the Hartree energy. Consequently, the exchange energy refers 
to the interaction between an exchange hole and the electron density over a certain range. Electrons with 
different spins can share the same orbit, but their mutual negative charge forces them to stay at a distance 
from each other. This electronic adjustment results in reduced electron density around each electron, 
resulting in lower attraction energy. This phenomenon is called the correlation hole [12]. The XC hole 
is predominantly determined by the exchange component at higher electron densities because of the 
Pauli exclusion principle, which is more pronounced when electrons are close to each other. However, 
at lower electron densities, the correlation component becomes more important and becomes 
comparable to the exchange component. Since a significant portion of the kinetic energy and the long- 
range Hartree energy are treated separately, the remaining XC energy can be reasonably approximated 
as a local or semi-local function of electron density. Furthermore, the shape of the XC hole is assumed 
to be spherically symmetric in three dimensions. The local XC energy per electron is therefore the 
electrostatic interaction energy between an electron at position 7 and the XC hole density at position Tr’ 
expressed as: 


—==> 


dr (eq.26) 


1 —- — Ny (FF) 
E,,[n] = si n(r)dr f [F-7'| 
Where n,,(7, 7’) denotes the average coupling factor of the exchange correlation hole, 


Ny(F,F') = fo nd (F,F')da (eq.27) 


Here A is the range —separation parameter. 
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Consequently, it becomes possible to define the exchange correlation density as follows: 
éxc[n@)] => f= vac r : dr" (eq-28) 


The exchange-correlation hole can be divided into two parts: the exchange hole (Fermi hole) and 


the correlation hole (Coulomb hole), 
Nye(7,7') =n, (7,7) + n,(7,7" (eq.29) 


The exchange hole, ny can be characterized using the Hartree-Fock expression for energy, 


E,[n@)] =5fn@)ar pee ar (eq.30) 


\7-7"| 


Therefore, the exchange-correlation functional can be formulated as follows: 


Ex-(n)] = fn@)excin@)]dr (eq.31) 


By clearly understanding the specific structure of the exchange-correlation density, we can approach 
the exact ground state. In the next section, we will examine various approximations aimed at estimating 


the effects of exchange and correlation. 
2.2.1.7.2 Exchange-Correlation Energy Approximation 


In order to determine the exact ground state of a system, knowledge of the exact exchange 
correlation functional is crucial. Unfortunately, obtaining the exact formula is a difficult task. To fill this 
gap, various approximations have been developed, including the local density approximation (LDA), 
the generalized gradient approximation (GGA), the hybrid GGA, the meta-GGA, the Tran-Blaha 
modified Becke-Johnson function (TB- mBJ) and more. In the following subsection, we will discuss 


LDA, GGA, TB-mBJ and HSE approximations in detail. 
a) Local Density Approximation (LDA) 


The local density approximation (LDA) [13] is considered one of the oldest, simplest and yet most 
important exchange-correlation functionals. In this approach, the actual exchange-correlation energy of 
a system is approximated by the exchange-correlation energy of a uniform electron gas with the same 
density. This uniform gas represents the only system for which the exact form of the exchange- 
correlation energy is known. The LDA relies solely on local density, resulting in a total energy 


expression: 


ELPAIn@)] = fn@ehe™[n@)]dr (eq.32) 
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In this context, ¢2™ denotes the exchange-correlation energy density of a homogeneous electron gas 
characterized by the electron density n(1*), while the exchange-correlation potential is obtained as an 
energy functional derivative, as shown below: 
ay) — SExe4In@) 
Vie [n@)] = eee (eq.33) 
LDA has been extended to include spin-polarized systems through the local spin-polarized density 


approximation (LSDA), which is expressed as follows: 
EXP Aly (7), nu (F)] = f n@yeke™[n, (7), nF) dr (eq.34) 


As already mentioned, the energy associated with the exchange hole n, can be obtained from the 


exact analytical formula of Dirac [14]: 


1 _ 
E,[n@)] = —2(2p34 ~ 28% (eq.35) 


4 412° r, Ts 
Where r, is the Seitz radius. 


The local density approximation (LDA) provides a reasonably accurate estimate of the exchange 
energy and typically remains within a 10% margin of error. However, the correlation energy, which is 
usually much smaller, is often overestimated by a factor of two. Interestingly, these two errors partially 
compensate for each other. LDA is found to be quite reliable in predicting ionization energies of atoms, 
dissociation energies of molecules, and cohesion energies, with accuracy ranging between 10% and 
20%. Surprisingly, LDA achieves an astonishingly precise accuracy of about 2% for bond lengths in 
both molecules and solids. Nevertheless, there are cases where LDA is inadequate, for example in 
systems that are strongly influenced by electron-electron interactions, such as heavy fermions. To 
overcome the limitations of LDA, the next natural progression is to include not only the density n(7) 
at a given point r but also the charge density gradient Vn(f) to account for inhomogeneous nature of 
true electron density. This advance leads to the development of the generalized gradient approximation 
(GGA) for exchange-correlation energy, a more sophisticated approach that increases the accuracy of 


the calculations. 
b) Generalized Gradient Approximation (GGA) 


Hohenberg and Kohn proposed an improvement to LDA [15, 16] by incorporating a higher-order 
density gradient expansion term, called the gradient expansion approximation (GEA). However, the 
GEA is unable to reproduce the essential features of the exchange correlation hole associated with the 


LDA and loses its physical meaning. This is evident in the violation of the sum rules and the lack of 
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restrictions on the negative exchange gap [17]. Despite these shortcomings, it provides a way to 
construct the exchange correlation hole for GGA [15] by using a real space cutoff of the GEA exchange 
hole. By introducing an analytical function called gain function F,, [n(7),Vn(1)], one can modify the 


LDA energy density to obtain the GGA exchange-correlation energy: 
ESS4[n(¥)] = [n@cSS4[n@]F x. [n@), Vn@)]d7r (eq.36) 


When the energy functional is expressed in integral form, the general expression for 


E&¢A[n(¥)] takes the following form: 


VEE In] = SE = nF) SEO! + een) (eq.37) 


If we include the consideration of spin variables, the expression for the energy functional in the 


framework of the generalized gradient approximation (GGA) becomes: 


ES in, ), mu) = f n@)e$E4 In, F), 1), Vn, *), Vny(F)] dr (eq.38) 


Perdew and Wang introduced a diverse set of GGA methods, one of which stood out and was later 
simplified to become known as the PBE functional, which bears the names of its creators: Perdew, 
Burke, and Ernzerhof [18]. GGA methods represented a significant advance over LDA because they 
took into account not only the local electron density but also its gradients. PBE in particular gained 
widespread recognition due to its simplicity and impressive accuracy. Its main features include 
considerations of the local electron density, its gradient and second-order gradients within the gain 
factors F(x) and F(c). This multifaceted approach has made PBE the first choice for researchers in 
various fields due to its ability to accurately describe a wide range of materials and their properties. 
However, it is important to recognize the limitations of LDA and GGA. Despite their advances, these 
methods are often inadequate when attempting to predict the band gaps in semiconductors or to capture 
the complicated electronic structures of highly correlated systems such as Mott insulators. In the field 
of DFT, the exchange-correlation energy plays an indispensable role. It encompasses the complex 
quantum mechanical effects arising from electron-electron interactions, an aspect not explicitly 
considered in LDA and GGA. Therefore, accurate estimation of this energy component is of paramount 


importance for ensuring the reliability of DFT calculations. 


In addition, a revised version of the PBE function, known as revPBE [19], emerged to address 
specific cases where improvements were required. It simplifies the mathematical formulation of the 
exchange energy while maintaining the accuracy of the original PBE function. Recently, a new GGA 
approximation for the exchange energy functional was introduced by Wu and Cohen aptly called it the 


WC approximation [20]. This innovation has made significant progress in improving the calculation of 
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structural properties in solid materials. What sets WC apart is its commendable computational efficiency 
and its unique property of being parameter-free, which represents a significant advance in the field of 
DFT. These developments highlight ongoing efforts to refine exchange-correlation functional in DFT, 
with a focus on improving its accuracy and applicability in the study of a variety of materials and 


phenomena in the fields of condensed matter physics and materials science. 


c) Modified Becke and Johnson Potential (mBJ) 


While the two previously detailed approximations provide accurate predictions for structural 
properties, their performance deteriorates when it comes to electronic properties, particularly bandgap 
widths. The gaps calculated using LDA or GGA are often significantly underestimated. To resolve this 
discrepancy, an approach was developed in 2006 by Becke and Johnson that aims to correct the energy 
gap values calculated by LDA and GGA [21]. A later development by Tran and Blaha in 2007 led to the 
modified Tran-Blaha potential (TB-mBJ), an adaptation of the Becke-Johnson exchange potential [22]. 


This change has the form: 


Ve) @) = eVBR@) + (Be -2)+ [5 vee (eq.39) 
Where, 
Ng(F) = >, |Wi\* is the electron density. 
t,(r) = ym |W; Vy; |? is the kinetic energy density. 
And VER) = — mG (1 — en %o) — 5 Xq(Fe*eM) (eq.40) 


is the Becke-Roussel potential [23], which was proposed to model the Coulomb potential arising from 


exchange holes. 


The term in Eq.40 is determined using the following equation: 


eee =i ES ut 

Where 
Q,(*) = =(V2n,(F) — 2yD,(@)) (eq.41) 

With 
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> =,  1|Vn,|? 
DoF) = 2tg() — re 


(eq.42) 


The y parameter is empirically determined by the least squares fitting of the Hartree-Fock 


exchange energy. 


The b, term is expressed by the following relationship: 


x3(F)e*o ; 


b,(r) = ( mL (eq.43) 


Or o represents the notation for spin. 


The change is reflected in the inclusion of the parameter “c” in the functional formula. It is worth 
noting that setting c = / restores the Becke-Johnson function. The proposed representation for “c” has 


the following form: 


a 10 Vn(7) > 1/2 
c= a+ B( [ar (eq.44) 


a and B are parameters that can be adjusted (with the values a = -0.012 and B = 1.023 (Bohr)!” in the 


WIEN?2K code [24]), while V.¢;; denotes the volume of the unit cell within the studied system. 


It is important to highlight that the mBJ exchange potential involves the exchange of electron holes. 
This exchange potential has been self-consistently integrated into the WIEN2K code [24]. In addition, 
the correlation potential is given, which is determined using one of the GGA versions. In our calculations 


we used the PBE-GGA correlation potential. 
d) Heyd-Scuseria-Ernzerhof Functional (HSE) 


Hybrid functionals belong to a class of approximations for the exchange-correlation energy 
functional within the DFT. They combine a portion of the exact exchange from Hartree-Fock theory 
with the remaining exchange-correlation energy from other sources, which can be either an “ab-initio” 
or empirical. The exact exchange energy functional is expressed in terms of Kohn-Sham orbitals rather 
than density and is therefore given the name implicit density functional. Axel Becke introduced the 
hybrid approach in 1993 [25]. Incorporating Hartree-Fock (HF) exchange, also known as exact 
exchange, provides a straightforward method to improve the accuracy of various molecular property 
calculations. These properties include atomization energies, bond lengths, and vibrational frequencies, 
which are often inadequately, represented using basic “ab-initio” functionals [26]. A widely used variant 


is the Heyd-Scuseria-Ernzerhof function (HSE) [27]. 
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The Heyd-Scuseria-Ermzerhof (HSE) exchange-correlation function [27] uses a failure function- 
tested Coulomb potential to calculate the exchange energy component with the aim of improving the 
computational efficiency. The HSE functionals are known for their ability to improve the precision of 
traditional semi-local functionals such as PBE, particularly when dealing with semiconductor bandgaps. 
In addition, they offer lower computational cost compared to other hybrid functions. The exchange- 


correlation energy has the following formula: 
Ege BEM = Bx *S (wn) + (1 — aE E(w) + Byes (w) + ECP* (w) (eq.45) 


In this context, “a” represents the mixing parameter, while “w” serves as an adjustable parameter that 
determines the short-range nature of the interaction. Commonly used values of “a@ = 1/4” and “w = 0.2” 


(often referred to as HSE06) have shown favorable results for most systems. The HSE exchange- 


correlation functional simplifies to the PBEO hybrid functional when “w = 0”. Here "E eek (w)" stands 


for the short-range Hartree-Fock exact exchange functional, "EF Eee SEN and "E eee denotes the short- 


range or long-range components of the PBE exchange function and "E ie : (w)" represents the PBE [18] 


correlation function. 
2.2.2 Augmented Plane Wave (APW) Basis Set 


While the value of pseudopotentials in computing electronic structures is undeniable, there is a 
growing trend to explore alternative methods for certain specific reasons. While pseudopotentials are 
highly efficient, they have limitations in capturing information that is inherently encoded near the atomic 
nucleus. In this context, the concept of Augmented Plane Waves (APW) is gaining importance and is 
proving to be a promising alternative. In the coming sections, we begin by introducing the Augmented 
Plane Wave (APW) method. We will then examine its improvements, including the Linearized 
Augmented Plane Wave (LAPW) method, Augmented Plane Waves coupled with local orbitals 
(APW+LO), and the Full Potential Linearized Augmented Plane Wave (FP-LAPW) approach. 


2.2.2.1 Augmented Plane Wave (APW) Method 


The Augmented Plane Wave (APW) method, originally developed by Slater [28] in 1937 to solve 
the Kohn-Sham equations, is based on a fundamental concept: the division of space into two different 
regions, as shown in Figure 3. In the inner region we encounter non-overlapping spheres, each centered 
on atomic locations, affectionately known as “muffin tin” spheres (abbreviated as MT spheres). These 
spheres have a characteristic radius called Rm. Within these limits, electrons assume a relatively free 
state that we can describe with plane waves. Here the potential exhibits spherical symmetry and the 
wave functions naturally assume radial shapes. Conversely, in the second region we are faced with what 


we call the “intermediate region’, which represents the unoccupied space woven between the muffin tin 


80 


Chapter 2 Theory and Computational Techniques 


balls. Within this region, electrons form strong bonds with the atomic nuclei, which resemble the 
behavior of free atoms. The potential remains constant and wave functions take the form of plane waves 
[29]. This fascinating duality in treatment allows the APW method to effectively capture the complex 


interplay between the behavior of free and bound electrons within a crystalline structure. 


Se ee a ee 


Figure 3. Illustration of muffin tin spheres and the interstitial region where potential singularities occur 
at the core positions. The Augmented Plane Waves (APW), shown in red lines, must transition smoothly 
exactly at the boundaries of the spheres and coincide with the plane waves in the interstitial region, 


shown by the blue line. This image was taken from the FLEUR user manual [30]. 


According to the Muffin-Tin approximation, the wave function wr) takes the following specific 


form: 


wr) = Ym Ami MY in®) r< Rint 
wi = aD ,€6 ei(G+KF > Rin (eq.46) 


Where 2 is the volume of the unit cell, Aj, and Cg are development coefficients, Yj_(1) is spherical 


harmonics and U{(f) is the regular solution of the radial part of the Schrédinger equation, given by: 


d? 1(1+1) > — 
{- ++ - £}r U,@ =0 (eq.47) 

In this equation, V(r) denotes the spherical component of the potential enclosed in the sphere, often 
referred to as the muffin tin potential. Meanwhile, E, represents the energy of linearization. These radial 
functions defined by this equation have orthogonality with respect to each nuclear eigenstate within the 
sphere. It is worth noting that this orthogonality property disappears at the spherical boundary [31], as 


explained by the following Schrédinger equation: 
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d77U, d77U2 
dt2 1 ar 


(Ez — E,)rU,U2 = Uz (eq.48) 


Where U, and U2 are radial solutions corresponding to energies Fi and E2, respectively. 


Slater's method introduces a key concept: plane waves serve as solutions to the Schrédinger 
equation when the potential remains constant, as it does in the interstitial region. In contrast, radial 
functions emerge as solutions when dealing with a spherical potential, and this is also true when E, 
represents an eigenvalue. While this approximation proves to be extremely effective for materials with 
cubic or face-centered structures, its applicability for asymmetric materials is limited. To maintain the 
continuity of the function y(F) at the surface of the MT sphere, the coefficients Aj,, must be expressed 
as Cg coefficients associated with plane waves within the interstitial region. These coefficients can be 


expressed explicitly as follows: 
4 il = > # = => 
Aim = RT aR CoJi(K + G|Rimt)Yim(K + G) (eq.49) 


The reference point is set at the center of a sphere with a radius of Rur. Consequently, the 
coefficients are expressed as plane wave coefficients, denoted Cg. The energy parameters Cg are 
referred to as coefficients of variation in the APW method. This alignment ensures that the functions 
indexed by G are compatible with the radial functions within the spheres, leading to the concept of 
Augmented Plane Waves (APW). These APW functions specifically serve as solutions to the 
Schrédinger equation for energies represented by E, alone. Therefore, it is imperative that the energy 
value E; matches that of the G-indexed energy band. This condition means that energy bands (for a 
given k-point) cannot be obtained by a simple diagonalization process. Instead, it becomes necessary to 


treat the secular determinant as a function of energy to decipher the complicated energy band structure. 


The inherent challenge within the APW method arises from the presence of the parameter E7 in the 
denominator of the U; function (Rut). This parameter can take the value zero at the surface of the MT 
sphere, resulting in a discontinuity between the radial functions and the plane wave functions. To address 
this problem, several modifications to the APW method have been introduced, the most famous being 
developed by Koelling [32] and Andersen [33]. Andersen's modification involves a novel approach. The 
aim is to represent the wave functions within the spheres by creating a linear combination of the radial 


functions U;(1) and their derivatives with respect to the energy Ey. This representation is formulated as 
a dU 

follows: U,(r) = ore which gave rise to the LAPW (Linearized Augmented Plane Wave) method. 
l 


This adjustment effectively bridges the gap and ensures a seamless transition between radial and plane 


wave functions, improving the overall accuracy and versatility of the method. 
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2.2.2.2 Linearized Augmented Plane Wave (LAPW) Method 


Once the fundamental functions and their derivatives are continuous and adapt to their status as 
radial functions with a fixed energy E;, we introduce a more adaptable and precise framework for 
understanding the band structure of solids - the LAPW (Linear Augmented Plane Wave) method. The 
LAPW method was a remarkable success, especially after Andersen's refinements of linearization 
techniques. Below we outline some key principles of the LAPW method. More details can be found in 
the work of Singh [34]. Within the LAPW method, especially in the atomic MT region, wave functions 
take the form of a linear combination involving radial solutions U;(7)Y;(7) and their derivatives, 
denoted as U;(7)Y;m(#), related to energy fluctuations. These functions U are defined according to Eq. 


48 from the APW method. In addition, it is important that the derivative functions U satisfy a certain 


condition: 


ad? Ut = po 7 
{- a2 + _ ) + Vcr) 7 E,|r U(r) = TU, (*) (eq.50) 

LAPWs are plane waves in the intermediate zone of the unit cell that reach the numerical radial 
functions inside the spheres, under the condition that the functions and their derivatives are continuous 


at the boundary as described in Eq.51: 


Yim|[AmU i) a Bim Ui) |Yim®) r < Rt 
T) = ais = 51 
yr) mc 7 >Rmt (eq.51) 


Where the coefficients By, corresponding to the function U;(7) are of the same type as the Ayn 


coefficients. 


As with the APW method, LAPW functions are plane waves only in the interstitial region. Within 
spheres, LAPW functions are more suitable than APW functions. If EF, differs slightly from the band 
energy E, a linear combination will actually reproduce the radial function better than APW functions 
consisting of a single radial function. Consequently, the U; function can be developed as a function of 


its derivative and the energy E: 
U,(E,r) = UE, Fr) + (E — E)u,(E,F) + OE — E,)* (eq.52) 
So that (E — E,)? represents the energy squared error. 


While this approach in the LAPW method guarantees the continuity of the wave function at the MT 


sphere boundary, it is worth noting that these calculations may experience a reduction in precision in 
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terms of computational accuracy compared to the APW method. The LAPW method leads to an error 
in the wave functions on the order of (E-E))? and another error in the band energies on the order of (E- 
E))*. Despite these errors, the LAPW wave functions serve as a robust basis for capturing all valence 
bands within a relatively wide energy range, all with a single energy value E;. In cases where this is not 
possible, splitting the energy window into two parts provides a significant simplification compared to 
the APW method. When U; takes a zero value at the surface of the MT sphere, its derivative U, is 
generally non-zero, which effectively solves the continuity problem at the surface of the MT sphere 
within the LAPW method. Takeda and Kubler [35] introduced a generalized version of the LAPW 
method using N radial functions and their (N-/) derivatives. Each radial function is associated with its 
unique parameter Ey, effectively mitigating the error associated with linearization. The default LAPW 
method is called when N=2 and El, is close to Elz. However, for N>2 these errors can be further reduced. 
Unfortunately, using higher order derivatives to achieve convergence requires significantly longer 
computation times. Singh [36] modified this approach by introducing local orbitals into the basis set 
without increasing the plane wave energy limit, which provided a more efficient strategy to deal with 


convergence while maintaining computational efficiency. 
e Role of Linearization Energies E; 


The extended wave functions U; and U, must satisfy the condition of orthogonality with respect to 
the core states within the MT sphere. However, this condition is only met if there are no nuclear states 
with the same quantum number “7”. Consequently, there is a risk of fusion of semi-nuclear states with 
valence states. The problem of non-orthogonality between certain nuclear states is still not addressed by 
the APW method, and the introduction of the LAPW method requires careful selection of the energy Ei. 
Therefore, an adjustment of £; is required to facilitate the calculation. In such cases, the ideal solution 
is to apply local orbital expansion. However, it is important to note that this option may not be available 
in all software programs. In such cases, selecting the largest possible MT sphere radius is crucial. 
Ultimately, the different EH; values should be defined independently of each other. Energy bands 
correspond to different orbitals. For precise electronic structure calculations, E; should be chosen as 


close as possible to the energy of the band when both have the same quantum number “7”. 
2.2.2.3 Development of LAPW into Local Orbitals 


The primary goal of the LAPW method is to obtain precise band energies close to the linearization 
energies E; [33]. For most materials it is sufficient to choose these energies near the center of the band. 
However, this approach is not always feasible because there are materials for which a single E; value is 
not sufficient to calculate all energy bands. This is particularly relevant for materials with 4f orbitals 


[37, 38] and transition metals [39, 40]. It highlights the fundamental problem of the semi-nuclear state, 
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which lies between the valence and the nuclear state. To address this challenge, multiple energy 


windows are often used or, alternatively, local orbital expansion is used. 
2.2.2.3.1 LAPW+LO Method 


The development of the LAPW method to local orbitals involves modifying the orbitals within their 
basis set to eliminate the need for multiple energy windows. This is achieved by introducing a third 
category of basic functions. The underlying principle is to treat all energy bands within a single energy 
window. Singh [36] articulated local orbitals as a linear combination of two radial functions 
corresponding to two different energies, along with the derivative of one of these functions with respect 


to the energy: 


[AimUi(T, E11) + Bim Ui(F, E11) + ComU(F; E12) ¥in®) 7 < Rint 


eq.53 
0 T>Rmt (9:52) 


p(F) = 
The coefficients Cj, have the same nature as the previously defined coefficients Aj, and Bin. 


A local orbital is defined for a given / and m and also for a particular atom, since within the unit 
cell all atoms must be taken into account, not just the inequivalent ones. In addition to addressing the 
half-core states, local orbitals can also be used to improve the basis set for conduction bands. This 
improvement to the LAPW method has been crucial to its widespread success, as it extends the 


applicability of this linearization approach to a significantly broader range of compounds. 
2.2.2.3.2 APW+LO Method 


The problem encountered with the APW method was the energy dependence of the entire basis set. 
While the LAPW+LO method managed to eliminate this energy dependence, it came with the 
disadvantage of using a larger base, which imposed significant limitations on both the APW and 
LAPW+LO methods. Sjésted, Nordstrém, and Singh [41] achieved a significant improvement by 
developing a hybrid basis that combines the advantages of the APW and LAPW+LO methods. This 
approach is referred to as APW+LO and requires only a slightly higher plane wave cutoff energy 
compared to the APW method. A standard APW basis is used and U;(7) is taken into account for a fixed 
energy £7), thereby retaining the advantages of linearizing the eigenvalue problem. However, since a 
fixed energy basis does not provide a satisfactory description of the eigenfunctions, local orbitals are 


introduced to ensure variational flexibility in the radial basis functions. 
An APW+LO basis is defined by the combination of two types of wave functions: 


e Plane waves (APW) with a set of fixed energies Er: 
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Yim|[AmUi@) + Bim Ui) |Yim®) r < Rmt 
(2 or 2 (eq.54) 
. ) CoeiGtr r > Rint q 
2L4G 


e Local orbitals different from those of the LAPW+LO method, defined as: 


YimlAmUi a E}) aS Bim U(r, E)|Yim®) r < Rt 
a (eq.55) 


Tas 


p= 


A mixed basis can be used in calculations, including both LAPW and LAPW+LO components, 
taking into account different atoms and even different values of the quantum number '/'. Typically, this 
approach assigns orbitals that exhibit slower convergence with the number of plane waves, such as the 
3d-states in transition metals or atoms with a smaller MT sphere size to the APW+LO basis. Meanwhile, 


the remaining components are described on the LAPW basis [42]. 
2.2.2.4 FP-LAPW Method 


The Full Potential Linearized Augmented Plane Wave (FP-LAPW) method [43] does not make any 
approximations regarding the shape of the potential or the charge density. Instead, they are developed 
into lattice harmonics within individual atomic spheres and represented as a Fourier series in the 
interstitial regions. This property is the origin of the method's name, “full-potential”. In particular, this 
approach ensures the continuity of the potential at the surface of the MT sphere and can be expressed in 


the following comprehensive form: 


: Vim( T)¥ im) 7 <R 
V(®) = dtm = Yim aE (eq.56) 
Lk Vxe'* r>Rme 
Likewise, the charge density develops in the following way: 
: Nim (EY im (E) TF <R 
(3) = pee ) Im a mt (eq.57) 
DK Nx e'* r>Rme 


In addition to the basic principles of the APW method, several variants and extensions have been 
developed. These include the Soler-Williams formulation of the LAPW [44], ELAPW [45, 46] and 
QAPW [47, 48] methods, each of which provides improved accuracy, flexibility, precision and 
computational efficiency compared to the standard LAPW method. This makes them valuable tools for 


electronic structure calculations for various types of materials and systems. 


Furthermore, when dealing with low-dimensional systems, it is possible to extend the division of 


the unit cell to explicitly include semi-infinite vacuum regions, each with tailored plane wave 
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amplification [49]. This approach enables efficient calculations for low-dimensional systems such as 
surfaces and thin films. In addition, an extension to one-dimensional setups was developed for 
addressing atom chains [50]. A large number of software projects have been developed to implement 
the LAPW method and its many variants. Some notable examples of these software codes are: Elk [51], 
Exciting [52], Flair [53], FLEUR [30], HiLAPW [54] and WIEN2k [42]. These software packages 
provide valuable resources for researchers working on the calculation of electronic structures. They 
provide a comprehensive toolbox for exploring the properties and behaviors of various materials, 
increasing the versatility and effectiveness of computational investigations in this area. Software 
programs such as VASP [55] and Quantum Espresso [56] adopt the pseudopotential approach alongside 
plane wave and projector-augmented waves (PAW). This approach is characterized by its user-friendly 
properties, requires fewer user control parameters and provides a streamlined development process. The 
use of plane waves also contributes to the simplicity of the code by simplifying various expressions in 


the software. A basic description of this method follows. 


2.2.3 Projector Augmented Wave (PAW) Method 


Several approaches have emerged to solve the Schrddinger equation. Among the prominent 
electronic structure techniques, we can divide them into two main groups: (i) linear methods, 
exemplified by the previously discussed Augmented Plane Wave (APW) method by Andersen [28], and 
(ti) the pseudopotential method, which is based on norm-conserving pseudopotentials of the first 
principle, which were developed by Hamann, Schluter and Chiang [57]. The pseudopotential method, 
in combination with a plane wave basis set, offers the advantage of formal simplicity. However, when 
it comes to first-row elements or systems containing d or f electrons, even pseudopotentials can become 
very “hard”. In such cases, practicality requires the use of either extensive or complicated basis sets 
instead of plane waves. Furthermore, the treatment of semi-nuclear states as valence states, which is 
often required for early transition metal and alkaline earth metal elements, leads to the use of hard 


pseudopotentials, which limits their transferability. 


Vanderbilt's ultrasoft pseudopotentials [58, 59] significantly alleviated this problem by relaxing the 
norm conservation condition typically imposed on the pseudopotential approach. This approach also 


allows economical treatment of first-row elements and transition metal. 


Car and Parrinello introduced a combination of density functional theory and molecular dynamics 
techniques [60] in which both electronic structure and atomic dynamics are solved simultaneously using 
Newton's equations. This approach simplifies structure determination and facilitates the dynamic 
evolution of atomic structures. The Car-Parrinello method was originally used in the context of the 


pseudopotential for plane waves. However, there is growing interest in applying this technique to all- 
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electron (AE) methods, which enable efficient handling of first-row elements and transition metals. AE 
methods provide insights into the behavior of the wave function near the nucleus, a feature not provided 
by the pseudopotential approach and which is crucial for various experimental techniques. Nevertheless, 
it is important to emphasize that up to this point, no low-energy molecular dynamics simulation [61] 


that rivals the quality of simulations using the pseudopotential approach has been performed. 


In this context, the concept of “Projector Augmented Wave” (PAW) functions was developed. 
Bléchl introduced the “Projector Augmented Wave” (PAW) method [62], which combines and extends 
the advantageous aspects of pseudopotential and LAPW methods through the use of projector functions. 
Basically, the PAW method matches partial wave functions derived from isolated atoms with pseudo 
partial waves, using specially designed projector functions that closely mimic the all-electron solution. 
In this framework, PAW datasets are generated that allow the recovery of full core electron wave 
functions, similar to the LAPW approach but with significantly reduced computational effort, 
approaching the efficiency of pseudopotentials (the basic concept of the PAW method is schematically 
presented in Figure 4). In simpler terms, valence wave functions exhibit fast oscillations near ionic 
nuclei due to the need to maintain orthogonality to the core states. This presents a challenge because 
accurately describing these wave functions requires a large number of Fourier components or, in the 
case of grid-based methods, an exceptionally fine grid. The PAW method addresses this problem by 
converting these rapidly oscillating wave functions into smoother counterparts, making them more 
computationally tractable. In addition, it provides a way to calculate the properties of all electrons based 
on these smoothed wave functions. This approach is somewhat similar to a transition from the 


Schréddinger picture to the Heisenberg picture. 


(b) 
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Figure 4. A schematic representation of the PAW method is shown. (a) Pseudo-quantities are 
established using a uniform, flat wave grid covering the entire space. (b) Pseudowave functions are 
reconstructed inside spheres and the respective single-center terms are subtracted. (c) The all-electron 
wave functions are reconstructed and the corresponding single-center energies are included. After Ref 


[63]. 
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In the Hilbert space, the wave functions (AE) are transformed into a new, so-called pseudo (PS) 
Hilbert space. The linear transformation denoted by ¢ converts the hypothetical pseudo wavefunction 


yp) into the real all-electron wavefunction |): 


lp) = FI) (eq.58) 


The “all-electron” wave function, which represents a single-body Kohn-Sham wave function, 
should not be confused with the many-body wave function. To ensure that {y) and |y) differ mainly in 


the regions around the ion nuclei, we can express it as follows: 


C=1+ DRG, (eq.59) 


The non-zero @,. values of are restricted to certain spherical augmentation regions Q,. surrounding 
the atom R. It is advantageous to represent the pseudowave function around each atom as an expansion 


of pseudopartial waves: 
|) = Lil@,) ¢; within Q, (eq.60) 


Due to the linearity of the operator ¢, the coefficients c; can be expressed as an inner product with 


a set of projector functions, denoted as |p;): 
ci = (p;|p) (eq.61) 
Where (pild,) = by. 
The all-electron partial waves, denoted as |@;) = {|@,), are usually chosen as solutions of the 
Kohn-Sham-Schrédinger equation for an isolated atom. The transformation ¢ is thus characterized by 


three components: a set of all-electron partial waves |@;), a set of pseudo-partial waves |,), and a set 


of projector functions |p;), and so it can also be expressed explicitly as : 


€=1+ DY; (1hi) — |b:)) (il (eq.62) 


Using this transformation, we can derive the all-electron wave function (AE) from the pseudo wave 


function (PS) by: 


Ip) = |b) + Xi (1s) — |61)) (pil) (eq-63) 


Outside the augmentation regions, the pseudo partial waves agree with the all-electron partial 
waves. Within the spheres, they can take any smooth continuation, including linear combinations of 


polynomials or Bessel functions. 
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2.3 Computational Formulation 


The main goal of practical DFT formulations is the numerical solution of the Kohn-Sham equations 
for systems consisting of multiple atoms. Density functional theory (DFT) calculations have quickly 
become a standard tool for a variety of materials modeling applications in physics, chemistry, materials 
science, and various engineering disciplines. The development of open-source electronic structure code 
packages, mainly based on the LAPW method or the pseudopotential method, has significantly enriched 
the toolbox for conducting detailed ab-initio materials studies. These joint efforts have produced well- 
structured shared codes that incorporate numerous state-of-the-art methods. A crucial aspect of code 
development is the validation process, with most collaborative teams incorporating internal testing as 
an integral part of their development processes. Additionally, the availability of multiple independently 


developed codes provides opportunities for extensive testing and validation. 


In this dissertation, we leverage the capabilities of two commercial software packages, WIEN2k 
[42] and VASP [55], to perform our DFT calculations. This section provides a detailed introduction to 


these two codes. 


2.3.1 WIEN2k Code 


The WIEN2K simulation code comes from the Institute of Materials Chemistry at the Vienna 
University of Technology and was first developed in 1990 by P. Blaha, K. Schwartz., P. Sorintin and 
S.B. Trickey [24, 64]. Since its inception, this code has been continually revised and updated multiple 
times, resulting in various versions known by their year of release, such as WIEN93, WIEN95, WIEN97, 
etc. In our work, we used the 2014 WIEN2K version to study the 3D bulk structures. 


WIEND2K is a computational software for calculating the Schrédinger wave equation for periodic 
materials in three-dimensional (3D), two-dimensional (2D) and one-dimensional (1D) configurations. 
It is coded primarily in FORTRAN90 and runs on the LINUX operating system. It consists of a 
collection of standalone programs linked together via C-SHELL SCRIPT. These programs are 
responsible for performing electronic structure calculations in solid-state materials using the principles 


of density functional theory (DFT). WIEN2K supports a wide range of function types and uses the FP- 
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LAPW method. It performs self-consistent calculations for all electrons, including both core and valence 


electrons, resulting in exceptionally accurate results. 


The code's capabilities extend to fundamental electronic structure calculations, structure 
optimization methods, spin-polarized/spin-orbit-based systems, phonons, and the simulation of various 
spectroscopic methods. For applications such as X-ray absorption or electron energy loss spectroscopy, 
WIEN2K considers excitonic effects by introducing a core hole on the relevant atom, thus enabling the 
precise modeling of various spectral features. Furthermore, WIEN2K offers the opportunity to go 
beyond the limits of DFT by leveraging advanced many-body perturbation theories such as the GW 
approximation [65] and the Bethe-Salpeter approach (BSE) [66]. 


2.3.2 WASP Code 


VASP, or the Vienna Ab-initio Simulation Package, is a powerful computational tool written 
primarily in FORTRAN. It is designed to perform ab initio quantum mechanical molecular dynamics 
(MD) of essentially small atom systems (approximately up to 100-200 atoms) using Vanderbilt 
pseudopotentials [58, 59] or the PAW method [62] combined with a plane wave basis set. While the 
core methodology of atoms focuses on DFT, VASP also considers post-DFT corrections, including 
hybrid functionals (such as HSE [27], PBEO [26] or B3LYP [67]) that integrate DFT with Hartree—Fock 
(RF) exchange as well as advanced techniques such as the GW method for many-body perturbation 


theory [65] and dynamic electronic correlations within the random phase approximation (RPA) [68]. 


Originally, VASP was based on code developed by Mike Payne, which also laid the foundation for 
CASTEP [69]. Later, VASP was brought to the University of Vienna, Austria by Jiirgen Hafner in July 
1989. The main program was created by Jiirgen Furthmiiller, who came to the Institute for Materials 
Physics in January 1993 together with Georg Kresse. The further development of VASP is currently 
being driven forward by Georg Kresse. Recent extensions include the adaptation of methods commonly 
used in molecular quantum chemistry for use in periodic systems. The latest version of VASP is version 
6.4.2, which became available on July 20, 2023. In our work, we used VASP version 6.2.1 to study the 


2D monolayer structures. 
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A more detailed description of the basic programs of the WIEN2K and VASP codes, step 


calculations, and possible properties can be found in Appendix B. 
2.3.3 BoltzTraP2 Code 


BoltzTraP2, an acronym for Boltzmann Transport Properties, is a widely used software package 
known for its ability to simulate and calculate the thermoelectric properties of materials. This code uses 
a contemporary approach that leverages smoothed Fourier interpolation of periodic functions and 
Onsager transport coefficients for extended systems and enables the careful calculation of a material's 
semi-classical Boltzmann transport equation (discussed in Chapter 1). It relies on crucial inputs 
including band- and k-dependent quasiparticle energies, intraband optical matrix elements, and 
scattering rates. BoltzTraP2 interfaces seamlessly with WIEN2K and VASP and provides versatile 
access through both a command line interface (via the btp2 command line interface) and a Python 


module [70]. 


The approach is based on expressing the quasiparticle energies and their derivatives for each band 


as Fourier series representations: 


e Quasi-particle energies €, are represented as follows: 
Ex = Laan Urea €Xp(ik. R) (eq.64) 
Where A represents sets of symmetry-equivalent lattice vectors. 
e Derivatives of quasi-particle energies Ve, are expressed as: 
Vex = iXa Ca Urea R exp(ik. R) (eq.65) 


BoltzTraP was originally inspired by Shankland's concept [71], which proposed determining 
coefficients by minimizing a roughness function while ensuring accurate reproduction of the calculated 
quasiparticle energies. This requirement implies that the number of coefficients should exceed the 


number of calculated points. 


The derivatives can also be derived from optical intraband matrix elements [72]: 


Vex = —(WxlPlWx ) (eq.66) 
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In BoltzTraP2, the Shankland algorithm [71] is extended to ensure that the coefficients reproduce 
not only the quasiparticle energies but also their derivatives, as expressed in Eq.66. This involves 


minimizing the Lagrangian: 
1 z S 
T= 5XRCrPR + UelAnl Ek — Ex) + LadakValEx — Ex)] (eq.67) 


The Fourier coefficients (Cr) are adjusted according to the constraints and the Lagrange multipliers 
(A, and Ag,) are chosen to satisfy these constraints. Here, the index @ denotes the three Cartesian 
directions, indicating that each calculated derivative introduces three Lagrange multipliers, as in Eq. 66. 


Pr is defined by Pickett et al. [73] as roughness function as follows: 
R R 
Pr=(1-y rare are Ge (eq.68) 


BoltzTraP2 uses the rigid band approximation (RBA) to calculate transport coefficients. This 
assumption implies that changes in temperature or doping level do not alter the underlying band 
structure. Furthermore, BoltzTraP2 is often associated with the Constant Relaxation Time 
Approximation (CRTA), which states that the Seebeck coefficient and the Hall coefficient are unaffected 
by the scattering rate [74]. This allows their determination on an absolute scale, taking into account 
doping and temperature fluctuations within a single analysis. However, it is important to note that under 
CRTA only electrical conductivity (o) and electronic thermal conductivity (ke) are reported, both of 
which depend on relaxation time (t) as a parameter. In order to accurately estimate «l and derive precise 


values for o and ke independently of t, alternative methods are examined in Chapter 3. 


2.4 Conclusion 


In summary, density functional theory (DFT) is a powerful computational tool for solving electronic 
structures in many-body systems. In this chapter, we have introduced the theoretical foundations 
underlying density functional theory. This includes discussions of the various approximations 
commonly used to solve the Schrédinger equation, as well as descriptions of possible implementations. 
Additionally, we have provided a brief overview of various software codes designed to solve different 


types of problems. 
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Chapter 3 


Effect of Doping and CO-Doping on the SrS Bulk 


Properties 


The first part of this work has been published in Journal of Solid State Communications 361, 115060, 
(2023). The second part of this work has been published in Journal of Materials Science in 
Semiconductor Processing 165, 107684, (2023). 


n the first part of this chapter, we discussed the spin-resolved properties of pure bulk SrS (3D) 
doped with iron (Fe) at different concentrations: 0%, 12.5%, 25%, 50%, and 75% Fe. In the second 
part, we focused on the Fe concentration of 12.5%. Here, we conducted further investigations by 
co-doping the compound with three different alkali metals, namely Li, Na, and K. This allowed us to 
thoroughly investigate the optical and thermoelectric properties. We also compared the properties of the 


singly-doped and co-doped systems and evaluated the improvements achieved. 


3.1 Structural, Mechanical, Electronic, and Magnetic Properties of mono-doped SrS: Fe 
Alloys 


3.1.1 Introduction 


Against the backdrop of ever-increasing information needs and improved social connectivity in 
today's world, the 21“ century has witnessed an unprecedented rise in technological advancement. This 
progress has been particularly evident in the significant development of high-performance computing 
systems and mobile devices. These advances have facilitated the seamless processing of large data sets, 
providing a robust response to the ongoing need for comprehensive data collection and storage methods 


[1]. 


The discipline of spintronics, which utilizes both electronic charge and spin properties, is at the 
forefront of technological innovation. It has enormous potential to pave the way for a new era 
characterized by fast, logically controlled and energy efficient data storage solutions [2]. The fabrication 


of dilute magnetic semiconductors (DMSs) characterized by high Curie temperatures (Tc), thereby 
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enabling the control of magnetism at ambient conditions by carrier and gate voltages, represents a key 
milestone in the realization of an excellent candidate for semiconductor spintronic technology [3, 4]. 
Recently, Fe-based DMSs have attracted great scientific curiosity. They are at the forefront of research 
due to their tantalizing potential as ferromagnetic semiconductors, promising impressive properties such 
as elevated Curie temperatures, minimal power consumption, and suitability for high-speed spin-based 
devices [5]. This intriguing potential has led researchers to delve deeper into the intricacies of the 
ferromagnetic mechanism in Fe-based strain gauges. Their joint efforts reinforce the amazing promise 
of iron to not only achieve semimetallicity at room temperature, but also to exhibit significant magnetic 
moment and mechanical robustness - especially when integrated into the group of II-VI chalcogenides 


[6-12]. 


As an archetypal member of the II-VI chalcogenide family, strontium sulfide (SrS) has attracted the 
attention of both the scientific community and the industrial sector alike due to its properties, which are 
discussed in detail earlier in Chapter 1. Various theoretical and experimental aspects work has been 
carried out on SrS alloyed with different proportions of magnetic transition elements (with 3d-orbitals) 
(see Section 1.6 of Chapter 1), but none has focused on conducting a detailed study on Fe-doped SrS at 
different concentrations. Furthermore, systematic studies on their elastic and mechanical properties are 
lacking. In this section of the chapter, we used a first-principles approach to investigate how doping 
concentration affects the structural, mechanical, electronic, and magnetic properties of Fe-doped rock- 
salt SrS. Our testing covers concentrations of 0%, 12.5%, 25%, 50% and 75% to ensure thorough 


analysis. 


3.1.2 Computational Details 


The computations discussed in this section were carried out using the Full-Potential Linearized 
Augmented Plane Wave (FP-LAPW) method [13] and the WIEN2k code [14]. The Generalized Gradient 
Approximation (GGA) as proposed by Perdew, Burke, and Ernzerhof (PBE96) [15] was utilized to 
address the exchange and correlation effects in order to determine the structural and mechanical 
properties. In addition to GGA-PBE, we optimized the results for electronic and magnetic properties by 
applying the Tran-Blaha modified Becke-Johnson (TB-mBJ) method [16], which is renowned for its 


excellent alignment with experimental data. 


In the FP-LAPW method, two regions take up space. The atomic sites are surrounded by non- 
overlapping spheres with muffin-tin radii (Rr) in the first region, and the interstitial region (IR) lying 
between the spheres is the second region. Spherical harmonics are used to extend the basis functions, 
electron densities and potentials within the spheres up to a maximum angular momentum of /inax=10. 


They are expressed in a Fourier series with a Fourier charge density in the interstitial region expanded 
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to Gnax=16 (a.u.)~‘ and cutoff radius Rur.Kmax=8. The last parameter, Kinax, is the size of the largest 
wave vector used for the plane wave expansion of the eigenfunctions and controls the size of the basis 
set used. Rur stands for smallest muffin-tin radius. The studied Sr).xFexS alloys (x = 0, 0.125, 0.25, 0.50 
and 0.75) have muffin-tin radii for the constituent atoms that are 2.4, 2.3 and 2.1 atomic units for Sr, Fe 
and S, respectively. We used the following electronic configurations to represent the valence states of 
these atoms in our calculations: Sr [5s”], Fe [3d° 4s”], and S [3s 3p*]. Brillouin zone sampling for SrS 
and Sr)-xFexS, respectively, is performed using an (8 x 8 x 8) and (4 x 4 x 4) gamma-centered k-point 
mesh to ensure computational efficiency. In successive iterations, self-consistency is considered to be 


achieved when the integrated energy and charge difference is less than 10* Ry. 


The elastic, mechanical and anisotropic properties of the compounds were investigated using the the 
Charpin method in conjunction with the WIEN2K code [17]. Using this method, it is possible to obtain 
the independent elastic constants Ci1, C12 and C44, which are suitable for describing the elastic behavior 


of compounds with cubic symmetry and help in determining other mechanical parameters. 
3.1.3 Results and Discussion 
3.1.3.1 Structural Properties 


In an ab-initio calculation, determining the structural properties is a crucial step in gaining further 
insight into the microscopic properties of the material. This leads to the study of other physical properties 
such as electronic, magnetic, etc. To proceed with the calculation of structural properties for ternary 
alloys at different concentrations, we studied the properties of their bare SrS. This allows for a more in- 


depth comparison and better understanding of the alloy properties. 


As mentionned in Chapter 1, bare SrS crystallizes under ambient conditions in the cubic rock-salt 
structure (B1) with a space group Fm3m (Ne.225), in which the Sr atom occupies (0, 0, 0), and the S 
atom occupies (0.5, 0.5, 0.5) Wyckoff positions, respectively (see Figure 9(a) in Chapter 1). 


For an investigation that captures both low and high concentrations, we used the supercell approach 
to construct a computationally economical and time-efficient (2 x 2 x 2) face-centered cubic (fcc) type 
supercell containing 16 atoms. This resulted in a series of rock-salt compounds, Srj.xFexS (x = 0.125, 
0.25, 0.50 and 0.75), with one, two, four and six of the eight Sr atoms replaced by Fe atoms, which gave 
Sr7FeiSg, SreFe2Ss, SraFe4Sg and Sr2Fe6Ss, configurations, respectively. All four of these compounds 
retain a cubic structure and each have one of three possible space groups. In particular, Sr7FeiSg (x = 
0.125) and Sr4Fe4S¢ (x = 0.50) retain the fcc structure with space groups 225 (Fm3m) and 227 (Fd3m), 
respectively. On the other hand, both SreFe2Ss (x = 0.25) and Sr2Fe6Ss (x = 0.75) adopt the primitive 


cubic system with space group 221 (Pm3m). We considered two possible couplings between atoms: 
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ferromagnetic (FM), in which the spins align in the same direction, and antiferromagnetic (AFM), in 
which the spins are opposite. In the case of the AFM configuration, a larger body-centered cubic (bcc) 
supercell with 32 atoms (2 x 2 x 2) was used. Figure 1(a-d) provides a visual representation of the crystal 


structure for the singly-doped SrS: Fe systems with the stable ferromagnetic phase (FM). 


Figure 1. Supercell models showing the cubic structure of Srj..Fe,S compounds with different x 
values (a) x = 0.125, (b) x = 0.25, (c) x = 0.50 and (d) x = 0.75. The Sr atom is shown in gray color, 


iron in red color, and sulfur in yellow color. 


The procedure for determining the structural properties near equilibrium involves evaluating the 
total energy of the systems for different lattice parameter values using the empirical Birch-Murnaghan 
equation of state (EOS) [18] as described by: 


E(V) = Ey + BoVo haan ( 


4 
ei 1 Vv 1 


Bu) \ Vo By Vo Bod (eq-1) 


Where Epo, Bo, Bo, and Vo, respectively, stand for the minimum energy, bulk modulus, the first 


derivative of the bulk modulus, and volume at the equilibrium state. 


The bulk modulus is determined at the minimum of the E(V) curve by the following relationship: 


BpHayo— (eq.2) 
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The volume is determined as a function of pressure by: 


-1 
V(P) = Vo (4 4 ar \M (eq.3) 

During the optimization process, we performed calculations in both spin-unpolarized and spin- 
polarized states. The variation of the total energy with volume for the four ternary compounds in the 
three configurations—ferromagnetic (FM), anti-ferromagnetic (AFM), and non-magnetic (NM)—using 
the GGA-PBE is shown in Figure 2. Based on Figure 2, it is clear that the ferromagnetic phase is the 
most stable and therefore favorable for all compounds because it has a significantly lower total energy 


compared to the anti-ferromagnetic and non-magnetic phases. 
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Figure 2. Change in total energy with volume for Srj..FexS compounds (x = 0.125, 0.25, 0.50, and 


0.75) in spin-polarized and unpolarized states. 
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The results obtained from the E(V) curves - lattice parameters 'a' and bulk modulus 'B' of the 
equilibrium unit cell for our binary compound SrS as well as for the four ternary compounds in the FM 
stable phase are summarized in Table 1, which provides a comparison between our results and previous 
experimental and theoretical results for pristine SrS. It is worth noting that there are no results in the 
literature that can be compared with our results for the ternary compounds. Therefore, our results serve 
as a reference point for other theoretical investigations and highlight the need for experimental 


validation. 


Table 1 


The optimized lattice parameters a (A), bulk moduli B (GPa), and formation energies Ey 


(eV/atom) for binary SrS and ternary Sr,-xFexS (x = 0.125, 0.25, 0.50, and 0.75) compounds. 


Compound a B Es 

Our work Exp Others Ourwork Exp Others Our work Exp Others 
SrS 6.059 6.024 6.061° 48.273 58° 47.184° - - - 

6.050° 48.300° 

Sros7sFeo.12sS = 5.949 - - 50.244 - - -4,384 - - 
Sro.7s Feo.25S 5.852 - 5.7794 52.410 - 59.5224 -4,324 - 1.5584 
Sro.s0 Feo.s0S 5.612 - - 60.329 - - -4.301 - - 
Sro.25 Feo.75S 5.336 - - 73.687 - - -4.422 - - 


* Ref [19], ° Ref [20], ° Ref [21], ‘Ref [22]. 


Examination of the data from this table shows a slight deviation between the calculated lattice 
constant and bulk modulus of SrS compared to experimental values [19]. The PBE method shows a 
tendency to overestimate the cell parameter by about 0.5% and underestimate the bulk modulus by about 
18.3%. Such discrepancies are consistent with the well-known trend of the PBE method. However, our 
results show excellent agreement with theoretical literature values using the PBE exchange potential 
[20, 21]. Furthermore, the estimated values for the ternary compound Sro7sFeo2sS are largely in 


agreement with the theoretically proposed values by Hoat [22]. 


Replacing strontium with iron results in a gradual reduction in the lattice constant as the 
concentration increases from 12.5% to 75%. This suggests a reduction in the size of the structure due to 
the smaller ionic radius of the substituted dopant Fe?* (0.55 A) compared to the bare cation Sr** (1.26 
A). At the same time, there is an increase in the compression modulus, which means that the structures 
become harder as Fe concentration increases. This behavior is consistent with findings from a study of 


similar transition metal-doped SrS using the GGA-PBE function within the FP-LAPW approach [20]. 
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Thermodynamic process analysis allows us to assess the difficulties in doping ordered structures 
and obtain further information about the stability of the systems under consideration. Fe-doped SrS 
thermodynamic stability was evaluated by calculating the formation energy (Ey) with the following 


formula: 
E¢ = Edopea — XE? (Sr) — yE’ (Fe) — zE”(S) (eq.4) 


Egopea denotes the total energy of the doped systems. E (Sr), E’(Fe), and E(S) refer to the per-atom 


energies of bulk Sr, Fe, and S, respectively. The variables "x," "y," and "z" represent the amounts of Sr, 


Fe and S atoms present in the unit cell. 


The formation energy results are presented in Table 1, where a negative and lower value of 
Ey indicates a higher degree of stability and simplicity in the doping process. Notably, all compounds 
exhibit negative Ey values, confirming their thermodynamic stability and suggesting their feasibility for 
synthesis under ambient conditions. The negative sign means that when these compounds form, thermal 


energy is released into the environment, indicating an exothermic reaction. 


Regarding the Ef values, the order is as follows: Er (Sro.25sFeo.75S) < Er (Sto.g75Feo.1255) < Ey (Sro.75 
Feo.25S) <Ef (Sro.soFeo.soS), indicating that Sro.2sFeo7sS is the easiest to synthesize, with the system 
appearing most stable at x = 0.75. Comparing our Ey value for Sro.7sFeo.2sS (-4.324 eV) with the 
previously published value using GGA-WC (-1.558 eV) [22], we find that the estimated value using 
GGA-PBE is 21% lower. This discrepancy can be attributed to differences in both the XC energy 


function used and the computational inputs. 
3.1.3.2 Mechanical Properties 


After examining the structural properties, we looked intensively at the mechanical properties. This 
step allowed us to not only verify the mechanical stability of the compounds but also investigate the 
elastic properties in the polycrystalline state and study their anisotropic responses. In our dissertation, 
all of the compounds we examined showed cubic symmetry. This property allowed us to rely on only 
three independent elastic constants (C11, Ci2, and Cay) to thoroughly characterize their mechanical 
properties. Determining the elastic constant requires the establishment of three equations to be solved, 


which are generated by applying three different types of deformation: 


The first step involves calculating the bulk modulus B: 
1 
B=7(C11 + 2€y2) (eq.5) 
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The second step involves applying a tetragonal stress tensor at constant volume to derive Cui - C12: 


o 0 0 
06d ; 0 (eq.6) 
0 0 


(1-0)? 


E= 


Where € represents the strain and o represents the applied stress. 
Within this stress, the total energy is given by: 
E (a) = E (0) + 3(€y1 — Cy2)Vo0" + 0(0)3 (eq.7) 
Where: 


E (0): is the energy of the unstrained system. 


Vo: is the volume of the unstrained system. 


The third step involves applying a rhombohedral stress tensor at constant volume to calculate the 


C44 constant, which is given by the following expression 
- 111 
oS 111 (eq.8) 
111 


The elastic constants can finally be determined and must satisfy the convergence criteria of Born- 


Huang [23] to be in the mechanically stable states and are provided as follows: 


Cy1 — Cy2 = 0; C4, 20 


C44 + 2C42 = 0; C44 = 0 (eq.9) 


The obtained values for the elastic constants C11, Ci2 and C44 of our binary and ternary compounds 
with some theoretical comparison values for the binary compound are listed in Table 2. For the ternary 


compounds, the existing literature lacks experimental and theoretical studies for possible comparisons. 


According to the results reported in Table 2 and the stability criteria of Born [23], we find that the 
studied compounds met all the criteria for mechanical stability, indicating that they are mechanically 


stable. 
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Table 2 

The calculated values of elastic constants Cj (GPa), bulk modulus B (GPa), shear modulus G 
(GPa), Young’s modulus E (GPa), Pugh B/G ratio, Poisson’s ratio v, Zener anisotropy factor A, 
longitudinal v, (m/s), transverse v, (m/s), and average v,, (m/s) elastic sound velocities, Debye 


Op (K) and melting T,, (K) temperatures of SrS and Sr).xFexS (x = 0.125, 0.25, 0.50, and 0.75) 


compounds. 
Parameters SrS Sros7sFeoi2sS — Sro.7s Feo.25S Sro.s0 Feo.soS Sro.25 Feo.7sS 
Cu 110.817 93.703 95.651 112.295 145.268 
100.4° 
113.9! - - - - 
Cr 19.049 28.448 30.823 34.381 37.706 
17.28 
15.75° - - - - 
Cag 32.811 12.970 16.509 26.983 8.396 
26.88° 
30.3! - - - - 
B 49.638 50.201 52.432 60.352 73.671 
G 37.535 18.959 21.705 31.268 19.576 
B/G 1.322 2.647 2.415 1.930 3.763 
E 46.935 50.517 57.219 79.989 53.949 
Vv 0.198 0.332 0.318 0.279 0.377 
0.715 0.397 0.509 0.692 0.156 
VI 3281 3227 3326 3631 3466 
Vt 1240 1617 1718 2010 1535 
Vin 1576 1814 1923 2239 1733 
Op 251 228 246 298 243 
Ty 1198 1107 1118 1217 1412 


© Ref [24], ‘ Ref [25], Ref [26]. 


Notably, C11, which refers to unidirectional compression along the [100] direction, is higher than 
C44 which represents shear deformation in the [010] direction. This suggests that these compounds are 
more susceptible to deformation under shear deformation than under unidirectional compressive 
deformation. The trend observed in the C1; values shows a decrease in the transition from the binary SrS 
compound to the ternary alloy with 12.5% Fe concentration. There is then a steady increase in the Ci 
values as the Fe concentration increases from 12.5% to 75%. In contrast, Ci2 shows a steady increase 
with increasing concentration. On the other hand, when moving from SrS to 12.5% Fe, Cu, initially 
decreases and then experiences an increase up to 50% Fe before decreasing from 50% to 75% Fe 


concentration. This highlights the significant influence of Fe concentration on the elastic constants of 
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the single crystal and implies that a higher Fe concentration reduces the shear along [010], which could 


affect the mechanical stability of the compounds. 


Our analysis suggests that the non-monotonic fluctuation of Ci; and C4 is due to the transition of 
fcc to primitive cubic symmetry. Yang et al. [27] observed a similar behavior in C44 when they studied 
the influence of 3d-4d element concentration on the elastic properties of an Mg solid solution using first- 
principles calculations. However, they did not provide an explanation for the anomaly that accompanied 


the sudden increase in C44. 


Understanding the mechanical properties of a compound not only allows assessment of its 
mechanical stability, but also enables comprehensive quantification of properties related to mechanical 
behavior, such as: hardness, brittleness, ductility, stiffness, chemical bonding and anisotropy. In 
practical applications, materials are often in a polycrystalline state. However, ab-initio studies typically 
focus on materials in their single crystal form. Voigt, Reuss, and Hill [28-30] have proposed 
approximations to predict the elastic behavior of polycrystalline materials based on the elastic constants 


calculated for single-crystal counterparts. 


The hardness of a material is determined by its polycrystalline bulk moduli (B) and shear moduli 
(G). The compression modulus is estimated using Eq.5, while the shear modulus is defined by the 


following relationships using the Voigt-Reuss-Hill approximations [28-30]: 


Gy = cue (eq.10) 


— _5(€11-C12) C44 
fx= 4€44+3(Cy1—Cy2) fea) 
Here, Gy represents the Voigt shear modulus corresponding to the upper limit of the values of G, and 


Gp represents the Reuss shear modulus for cubic crystals corresponding to the lower values. 


The average of these two moduli represents the Hill shear modulus (Gy) and is given by the 
following equation: 


= Grt+Gy 


Gy 5 


(eq.12) 


The values of B range from 49.638 to 73.671 GPa and G range from 18.959 to 37.535 GPa, as listed 
in Table 2. This indicates that the pristine SrS with the lower B modulus value is the softest and that 
Sro.g75Feo.125S 1s less resistant to shear deformation compared to all other compounds. The trend for B is 
consistent and gradually increases with concentration. On the other hand, G can exhibit non-monotonic 
behavior, decreasing in the transition from pure to ternary alloy (from SrS to 12.5% Fe), increasing in 


the range of 12.5% to 50% Fe, and then decreasing sharply in the range of 50% up to 75% Fe. This 
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variation can be attributed to the non-monotonic elastic constants that contribute to such fluctuations. It 
is noteworthy that the B values calculated from the elastic constants agree closely with those obtained 


by the Birch-Murnaghan EOS fitting (see Table 1), confirming the accuracy of our calculations. 


The brittle/ductile behavior can be predicted by the Pugh ratio (B/G) [31], which relates the bulk 
and shear moduli to a critical value of 1.75. The material is considered brittle if the B/G value is less 
than 1.75; otherwise, it is ductile. As indicated in Table 2, the B/G ratio for pristine SrS is below the 
critical value of 1.75, classifying it as a brittle material, while for our four alloys it is above the critical 


value, classifying them as ductile materials. 


The stiffness of a material is quantified by the Young modulus (£), which establishes the 


relationship between the compression modulus and the shear modulus and is expressed as follows: 


__9BG 
~~ 36+B 


(eq.13) 


As Table 2 shows, the Young modulus actually exhibits a pattern similar to the shear modulus due 
to their interdependence. The addition of Fe to the SrS host matrix leads to an increase in stiffness, 
highlighting the reinforcing capacity of iron. A high modulus of elasticity indicates a material's ability 
to withstand stress. Of the compositions examined, the one with 50% iron content is expected to have 


the highest stiffness and the one with 12.5% iron content is expected to have the lowest stiffness. 


We also calculated Poisson's ratio (V), a quantitative measure of the type of chemical bond present 


in the material. This relationship can be expressed in terms of volume and shear modulus using Eq.14. 


_ 3B-2G 
~ 2(3B+G) 


(eq.14) 


If v exceeds 0.25, the chemical bond is said to be ionic; if it falls below 0.25, it is considered 
covalent [32]. According to the v values listed in Table 2, all ternary alloys have a v value in the range 
of 0.3. This implies that the bonding in the compounds Sro.g75Feo.125S, Sto.75Feo.25S, Sto.soFeo.s0S and Sro.25 
Feo.75S is mainly ionic. On the other hand, the value for binary SrS is below 0.25, indicating its covalent 


nature. 


The anisotropy of a material is determined using the Zener anisotropy factor (A). This factor gives 
a numerical indication of the extent to which a compound deviates from isotropy and is calculated using 


the following equation [33]: 


2044 
C11-Cy2 


(eq.15) 
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If the factor A is equal to one, this indicates perfect isotropy. When A is less than one, the material 
is harder in the [100] direction, and when A is greater than one, it is harder in the diagonal [111] direction 
[33]. As can be seen in Table 2, all compounds exhibit deviations from unity, showing their anisotropic 
nature with optimal stiffness in the diagonal direction [100]. The trend in A closely follows that of C44 
due to their interdependence. The greater the deviation from unity, the more pronounced the degree of 
anisotropy. This means that the iron-enriched compound (75% Fe) has the highest degree of anisotropy, 


while the 0% iron compound (SrS) has the lowest degree. 


To provide a more comprehensive understanding of the elastic anisotropy in the studied compounds, 
we used three-dimensional (3D) surface representations. These plots provide a visual representation of 
the variation of the polycrystalline modules across the main crystal directions. The specific mathematical 


formulations are listed below [34, 35]: 


1 
eo ait 2$42)(G+ + ) (eq.16) 
if: 
G = (S44 + 450J) (eq.17) 
5 = S11 — 2(S11 — Siz — 9) (GB+BG + GB) (eq.18) 


Where $44,542,544 are the elastic compliance constants; 1%, 1%, and [3 are the direction 


cosines; and J is quantified by: 
J = sin’0.cos?6 + 0.125. sin*@ (1 — cos 4@) (eq.19) 
6 and @ are Euler angles. 


Figures 3(a-c) and Figures 4(a-c) illustrate the 3D contour plots of the bulk, shear and Young's 
moduli for SrS and Sr1..FexS (x = 0.125, 0.25, 0.50 and 0.75), respectively. The obvious spherical shape 
in the 3D contour plots of linear compressibility of cubic systems, as shown in Figure 3(a) and Figures 
4(al-a4) confirm their isotropic nature. When examining the shear and Young moduli in Figures 3(b-c) 
and 4(b1-b4)/ (cl-c4), their 3D representations deviate from a perfect sphere, indicating anisotropic 
properties. This anisotropy is more pronounced in the compound with the highest iron content (75%). 
Notably, the Young modulus for all compounds peaks in the [100] direction and reaches a minimum in 
the [110] direction, which is consistent with the Zener anisotropy factor observations. However, it is 
important to note that the directional dependencies of the shear modulus are opposite to those of the 


Young's modulus. 
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Figure 3. The directional dependencies of three-dimensional (3D) contour plots of (a) bulk modulus 


(B), (b) shear modulus (G), and (c) Young's modulus (E) for SrS pristine compound. 


Figure 4. The directional dependencies of three-dimensional (3D) contour plots of (a) bulk modulus 
(B), (b) shear modulus (G), and (c) Young’s modulus (E) for Sr;..Fe,S compounds [(1) x=0.125, (2) 
x=0.25, (3) x=0.50 and (4) x=0.75]. 


Since mechanical behavior and thermal properties are highly correlated, the identified elastic moduli 


are mostly used to describe the melting point (T;,,.) and the Debye temperature (9p). 


The highest frequency of oscillation is related to an important physical parameter called the Debye 
temperature. From this temperature, important thermophysical parameters such as entropy, melting 
temperature and heat capacity can be calculated. According to the harmonic approach, the vibration 


frequency and the square root of hardness are directly related, and the Debye temperature can be used 
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to determine the hardness [36]. Using the data of the average speed of sound (¥,,,) and the elastic 


constant, the Debye temperature (8p) can be calculated in a conventional way [37]: 


oo= & EE ()P on 420 


Here, h/Kg represents the ratio of Planck's constant to Boltzmann's constant; n stands for the number 
of atoms in the cell; N4 is Avogadro's number; p denotes the mass density; M represents the molecular 


weight. Additionally, V,, encompasses both transverse (V;) and longitudinal (V,) velocities, outlined as 


[37]: 


Um = [5 (3+ 4) (eq.21) 


VI 


By applying the Navier equation [37] to the mass density, shear modulus, and bulk modulus, the 


transverse velocity (Vz) and longitudinal velocity (V,) can be obtained: 


5G 
v, = [Es = (eq.22) 


The heat sensitivity of a material is indicated by its melting temperature (T,,). Because they are 


stable over a wide range of temperatures, materials with a high melting point are highly desirable. To 
estimate the melting temperature of these materials, the empirical equation of Fine [38] is used, taking 


into account the elastic constant C)1. 
Tm(K) = 553K + (5.911K/GPa) C,, +300K (eq.23) 


The Debye temperatures and melting temperatures are shown in Table 2 and are graphically shown 
in Figure 5 for better visualization. It is obvious that the melting temperature has a nonlinear variation 
with the Fe concentration in the compounds and shows a similar trend with the elastic constant Ci, 
reflecting the direct correlation between them, as presented in Eq.23. The T,,, values obtained for all 
compounds are high, indicating their resistance to high-temperature treatment and their ability to 


maintain their solid-state structures at elevated temperatures. 


In contrast to the melting temperature, the Debye temperature shows a decrease with increasing 
concentration from 0% to 12.5%, followed by an increasing trend at 12.5% Fe and a subsequent 
decreasing trend at 50% Fe. This pattern is already evident in the elastic properties results, since the 
Debye temperature is an average of the sound velocity that correlates with the compression modulus 


and in particular with the shear modulus (see equations (20)—(22)). This suggests that the shear modulus 
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has a stronger influence on the variation of the Debye temperature. The compound Sro.soFeo.soS has the 


highest value, which indicates its relatively high intrinsic hardness. 


Unfortunately, direct comparison of our thermal and mechanical results with previous 
experimental work on this doped chalcogenide is not possible. However, it is interesting to note that our 


results agree with most previous first-principles calculations [39-41]. 


8p / T,, temperature (K) 


0 0.125 0.25 0.50 0.75 


Figure 5. The dependence of melting temperature (T,,,) and Debye temperature (9p) on Fe content for 


SrixFexS Compounds (x = 0, 0.125, 0.25, 0.50, and 0.75). 


3.1.3.3 Electronic Properties 
3.1.3.3.1 Electronic Band Structure Analysis 


After studying the effects of Fe doping on various properties and system stability, we will now 
consider how the introduction of Fe affects the electronic behavior of the SrS semiconductor. This is 
achieved through an analysis of spin-resolved band structures (BS), total density of states (TDOS) and 
partial density of states (PDOS). It is worth noting that previous studies have shown that bandgaps in 
materials are often slightly underestimated when using PBE-GGA, which is due to self-interaction errors 
[15]. An alternative is to use the mBJ exchange functional, which has been shown to provide reasonably 
accurate bandgap calculations that are close to experimental values. Figures 6 and 7 show the band 
structures along the high-symmetry points of the irreducible Brillouin zone for the host SrS and the 
ternary-doped Sr)..FexS compounds (x = 0.125, 0.25, 0.5, and 0.75, respectively) , calculated based on 


their equilibrium lattice constants. 
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Figure 7. Spin-resolved electronic band Structures of Srix.FexS compounds (x = 0.125, 0.25, 
0.50, and 0.75) with mBJ-PBE. 


According to Figure 6, SrS has a band structure typical of a semiconductor. In an mBJ calculation, 
the gap is indirect in the direction T-X with a value of 3.435 eV. Compared to other theoretical 
calculations, this value is closer to that of Hamidane et al. [20], with a difference of 0.135 eV, than that 
calculated by Kaneko and Koda [43], which has a deviation of 0.419 eV. The calculated gap value is 
underestimated compared to the experimentally measured value. It shows a difference of 26% compared 
to the experimental value of Sharma et al. [42]. However, this discrepancy arises from the use of 
different input parameters and the well-known tendency of DFT to underestimate energy gaps compared 
to experimental values. Our calculations compared to other theoretical and experimental results are 


summarized in Table 3. 


From Figure 7, which shows the spin-majority and spin-minority band structures for the ternary 
compounds, it can be seen that the introduction of Fe into the SrS host matrix results in a significant 
transformation of the electronic structure. In fact, the addition of 12.5%, 25% and 50% Fe content 
changes the nature of pristine SrS from a non-magnetic semiconductor to a ferromagnetic half-metal 


(HMF) and from ferromagnetic half-metal to magnetic metal at a higher degree of concentration (75% 
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Fe ). The semimetallicity in the first three concentrations arises from the fact that the majority-spin (left 
image) maintains the semiconducting behavior of SrS, with the Fermi level (Er) lying within the 
generated bandgap. Conversely, the minority-spin (right image) shows metallic behavior due to the 
overlap of some valence bands with the Ep, resulting in 100% spin-polarization. The Sro.g7sFeo.125S 
compound exhibits a direct bandgap with both the conduction band minima (CBM) and valence band 
maxima (VBM) located at the [ symmetry point with a value of 3.401 eV. However, the scenario 
changes drastically with increasing concentration, with the direct bandgap transforming into an indirect 
bandgap located at I'-X for Sro.7sFeo.25S, with a transition energy gap of 1.864 eV and at L-I’ for Sroso 
Feo.soS, with a transition energy gap of 3.539 eV. 


The energy gap (Ez) is defined as the energy difference between the maximum of the valence band 
(VBM) and the minimum of the conduction band (CBM). Thus, the E, is calculated for our compounds 
in the spin-up direction. The half-metallic gap (Gum) is another important parameter for characterizing 
half-metallic materials, which can be determined by the difference between the VBM and the Ef in the 
spin-up channel. The E, and Guo values for the three semimetallic ferromagnetic compounds are always 
given in Table 3. As indicated in Table 3, both E, and Gym values show nonlinear behavior with 
increasing iron content. The trend over the entire concentration range follows a specific order: a decrease 
from 12.5% to 25% Fe, followed by an increase from 25% to 50% Fe, and finally a complete 
disappearance between 50% and 75% Fe. This changing trend is consistent with observations made on 
previously studied TM-doped II-VI compounds [44, 45]. This behavior is attributed to the introduction 


of impurity states by Fe doping into the SrS bandgap, resulting in a sudden and significant narrowing. 


At a low Fe concentration (12.5%), the impurity states shift both the VBM and CBM slightly 
downward, resulting in a net bandgap reduction of 0.034 eV compared to pure SrS. This downward shift 
is attributed to the weakened interaction between the cation Sr and the anion S due to Fe doping. When 
the Fe concentration increases to 25%, the impurity states broaden, resulting in a significant downward 
shift in the CBM and a slight upward shift in the VBM. This results in a gradual reduction in the bandgap 
with a net reduction of 1.537 eV. The reduced bandgap value can be attributed to the combined effects 
of band modification due to the convergence of energy levels of the 3d bands of the Fe atom and the 3p 
bands of the S atom. This leads to an exchange interaction between these bands, as described in the 
following section, along with strain effects caused by the doping level. When the dopant concentration 
increases to 50%, the bandgap begins to increase gradually and reaches a maximum of 3.539 eV. The 
increased concentration of impurities leads to the formation of larger defects. Consequently, a higher 
defect concentration along with symmetry changes from fcc to primitive cubic may contribute to the 
observed broadening between the 3d-Fe and 3p-S bands (see the next section for more details). At 75% 
Fe, VBM and CBM follow a similar pattern to the band structure at 25% Fe, further reducing the 


bandgap until it finally disappears. In addition to the previously mentioned factors, the resulting 
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hybridization of 3d-Fe with 3p-S bands at the highest doping concentration could be a factor contributing 


to the sudden disappearance of the bandgap. 


Note that our calculated E, and Gym band gaps for Sro.7sFeo.25S show an underestimation compared 
to the previous theoretical values [22], with deviations estimated at 0.253 eV for E, and 0.571 eV for 
Guo. This underestimation is due to the differences in the approximation methods used (WC-based 
mBJ). An important point to consider is that both experimental [46] and theoretical [47] evidence have 
shown that a material's bandgap is directly correlated with its hardness. Consequently, the maximum 
bandgap observed in the Sro.s0 Feo.soS compound is consistent with the findings regarding its mechanical 


properties and Debye temperature. 


3.1.3.3.2 Total and Partial Densities of States Analysis 


To shed light on the role of atomic states in shaping semimetallic and metallic behavior and to 
further investigate the underlying reasons for the bandgap wobble, we performed comprehensive 
calculations of the TDOS and PDOS for pristine SrS and all Sri+.FexS compounds (x = 0.125, 0.25, 0.50 


and 0.75). The resulting curves for both spin configurations are shown in Figures 8 and 9, respectively. 


From the TDOS curves in Figure 8, it can be seen that there is no polarization between the electron 
densities of the states for up and down spin orientations of SrS. The density of states is symmetric and 
indicates the non-magnetic state of the compound, in contrast to the density of states for Sri..FexS, which 


is asymmetric (Figure 9). This is logically attributed to the doping effect with Fe. 
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Figure 8. The Spin-resolved total and partial density of states of rock-salt SrS with mBJ-PBE method. 
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Figure 9. Spin-resolved total and partial densities of states of Srj.Fe,S compounds (x = 0.125, 0.25, 
0.50, and 0.75) with mBJ-PBE. 


From the Sr curve (Figure 8), a significant bandgap is evident, indicating the semiconductor nature 
of this system. This is because the CB and VB do not intersect the Fermi level (Er), which is defined at 
0 eV. This bandgap arises mainly from the shift of the p-states of the S atom, which dominate the VBM, 


and the s-states of the Sr atom, which make a small contribution in both regions. 


Upon closer inspection of the TDOS curves of the ternary HMF compounds (Figure 9), a clear 
pattern emerges. In particular, in majority-spin, the 3d-Fe states and 3p-S states of Sro.g7sFeo.125S, 
Sro.7sFeo.25S, and SrosoFeo.soS, hybridize with each other, leading to the formation of a bandgap near the 
Er. This configuration gives these materials a semiconductor-like character. Conversely, in minority- 
spin, the 3d-Fe and 3p-S states intersect the Er, resulting in a distinct metallic character, confirming that 
these compounds are HMF with 100% spin-polarization at Er. In the Sro2sFeo7sS compound, both 
majority- and minority-spin channels exhibit metallic character that is due to the formation of 3d-Fe and 
3p-S states at the Er. This means that p-d hybridization is not particularly pronounced when the 


concentration of 3d-Fe electrons is relatively low. 


Analysis of the PDOS curves shows that for both Sro.g7sFeo.125S, Sto.7sFeo.2sS compounds, the VB in 
the energy range from about ~ — 4.14 eV to about ~ — 1.63 eV is mainly derived from 3d-Fe states, with 
a small contribution from 3p-S states. In the subsequent energy range from about ~-1.63 eV to about ~ 
— 0.42 eV, the contribution comes predominantly from 3p-S states. When the Fe content exceeds 25%, 
the contribution of 3p-S states becomes prominent alongside that of 3d-Fe states in the energy range 
from about ~ — 4.53 eV to about ~ — 0.74 eV for Sro.soFeo.soS and significantly ~ — 5.57 eV to about ~ 
— 0.05 eV for Sro.25Feo.75S. However, the involvement of 5s-Sr states is negligible for all compounds. 
Conversely, in the minority-spin channel, the lower part of the CB is mainly dominated by 3d-Fe states. 


These states start at 1.35, 1.02, 0.99, and 0.48-6 eV, respectively, and follow the order 0.125 — 0.25 > 
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0.50 — 0.75. States at the Er in the PDOS are predominantly attributed to the 3d-Fe states, with minimal 
contribution from the 3p-S states. These states play a crucial role in producing the ferromagnetic 
properties observed in the resulting compounds. The calculated PDOS values at the Er for the studied 
cases are summarized in Table 3. It is generally accepted that lower densities of states at the Er are 
associated with greater material stability [42]. This suggests that the compound with the highest iron 
content (75%) has the highest stability, which is consistent with the conclusions from the enthalpy of 


formation and melting point calculations. 


Table 3 


The calculated spin-majority bandgaps E, (eV), half-metallic bandgaps Gum (eV), and partial density of 
states (PDOS) at the Fermi level (states/eV) of SrixFexS compounds (x = 0.125, 0.25, 0.50, and 0.75). 


Compound Nature Location Eg Gum PDOS (Er) 
of gap Our work Exp Others Our work Exp Others 3d (Fe) 5s(Sr) 3p(S) 

SrS indirect [I-X] 3.435 4.328 + 3.300° - - - - - - 
3.016 - - - - - - 

Sro.s75 Feo.125S_ direct [T-T] 3.401 - - 0.431 - - 5.322 2.806x10* 0.005 

Sro7sFeo2sS indirect [X-I] 1.864 - 2.1174 0.090 - 0.6614 2.142 3.626x10* 0.008 

SrosoFeosoS indirect [L-T] 3.539 - - 1.059 - - 0.647 9.056x10* 0.086 

Sro.25 Feo.75S - - - - - - - - 0.432 1.848x10* 0.001 


" Ref [42], ° Ref [20], ' Ref [43], ‘Ref [22]. 


A more detailed understanding of the mechanism of interaction between the introduced Fe impurity 
states and the anion S states is provided by crystal field theory. The electron configuration of the Fe?* 
cation is [Ar] 3d°. The Fe?* ion forms six bonds with the ligand S atoms in an octahedral complex. 
According to Hund's rule and crystal field theory, each atom contributes two electrons, so a total of two 
electrons are added to the VB of the host semiconductor. The octahedral field created by the S anion 
after the introduction of Fe impurities splits the 3d-Fe states into two energy-different states: a triplet tog 
(dxy, dy and d,z) and a doublet eg (dz” and dx ay). Five of the six Fe electrons have majority-spin and 
occupy the states eg and tz, which are below the Er. In contrast, at the lower end of the conduction band 
the unoccupied eg states are in minority spin, and third of the tz, minority states are filled and raised to 
the Ep. Figure 10 shows a schematic diagram illustrating how the Fe** ion splits into its crystal fields in 


the octahedral environment. 
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Figure 10. Schematic diagram of the distribution of electrons in the eg and tzg levels of the Fe** ion 


situated in an octahedral crystal field. After Ref [50]. 


The main cause of the partially filled Fe-t2, states observed at the Er is the significant p-d 
hybridization. This hybridization also explains the two symmetrical states separation in the majority 
spin, which eventually causes ferromagnetism in the compounds under study to emerge and stabilize. 
The proposal of Zunger and Katayama [48, 49] on the ferromagnetism of octahedrally doped 


semiconductor systems with transition metals is similar to our observations. 


Examination of the Fe:3d direct exchange splitting energy Ax(d), determined by the difference (Ad| 
- Ad*t) between the respective up and down peaks, provides additional evidence for the existence of 
ferromagnetic alignment. This interaction arises from overlapping electron wavefunctions of 
neighboring atoms, resulting in an energy reduction when the electron spins align in parallel. This 
reduction in energy is called the direct exchange splitting energy. The calculated values are listed in 
Table 4. The high Ax(d) values confirm the consolidation of ferromagnetic alignment across all 


compounds. 


The spin-polarization P at the Er was determined and expressed by Eq.39 in Chapter 1. Table 4 
presents the results, which strongly validate the HMF nature of the compounds Srj_xFe,xS (x = 0.125, 
0.25, and 0.50), showing a magnetic spin-polarization of 100%. On the other hand, the compound Sri- 
xFexS (x = 0.75) exhibits a 33% polarization. This is obtained from 8% of the majority-spin states and 


4% of the minority-spin states that are both present at the Ep, strongly indicating that it is metallic. 
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Table 4 


Calculated values of exchange-splitting A,.(d), majority-spin N;(E;) and minority-spin N, (Ef) states 
at Ep, and the spin polarization P of Sri..Fe,S (x=0.125, 0.25, 0.50, and 0.75) compounds. 


Compound A,(d) N; (EF) N, (Ef) P 
(eV) (states/eV) (states/eV) 

Sro.875 Feo.125S 4.829 0 0.06 1 

Sro.7s Feo.25S 4.819 0 0.05 1 

Sro.so Feo.soS 6.748 0 1.33 1 

Sro.2s Feo.7sS 5.726 0.08 0.04 0.33 


3.1.3.4 Magnetic Properties 
3.1.3.4.1 Magnetic Moment 


In this section, we calculated the total magnetic moment of Sr).xFe,S at different concentrations (x), 
as well as the local magnetic moments of Sr atom, S atom, Fe impurity and the interstitial magnetic 
moment to clarify the cause of the induced magnetism in pristine SrS. Table 5 rearranges the collected 


data. 


The compounds Sro.37sFeo.125S, Sro.7sFeo.2sS and Sro.soFeosoS have been proven to be half-metallic 
ferromagnets because their total magnetic moment has an integer value of 4 us per Fe atom. The unfilled 
3d-Fe states are responsible for this measured magnetic moment. In this case, S has 3s” 3p* with two 
missing electrons in the outermost 3p shell, while Fe has an electronic valence configuration of 3d° 4s”. 
Two of the eight valence electrons are added to the dangling bonds of the anion (S) when incorporating 
Fe. The s of Fe electrons make no contribution to the magnetic moment because they are completely 
localized. The magnetic state is established by the remaining four d electrons of Fe, which are 
delocalized and fill the majority-spin while the minority-spin remains empty. The magnetic moment of 
the iron-rich compound (75% Fe) is characterized by a slight deviation from the expected integer value 


(48), a property that often occurs in metals. 


Small permanent local magnetic moments, which show a positive sign and are indicative of parallel 
spin arrangements (FM), were produced on the non-magnetic sites of the Sr and S atoms as a result of 
pd-type electron transfer. Additionally, Table 5 indicates that the Fe atom has the largest values and 
contributes the most to the compounds’ total magnetic moment, which decreases as concentration 


increases. The unfilled 3d electronic states of the magnetic atom Fe are responsible for this magnetic 
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moment. This pattern is consistent with previous research on SrS doped with transition metals at 
different concentrations [44, 20]. Furthermore, the interstitial sites positively contribute to the supercell's 
overall magnetic moment. This suggests that ferromagnetic interaction occurs at these sites between the 
Fe-3d states and the corresponding 5s-3p states of Sr and Sr. When comparing the predicted magnetic 
moments of Sro.7sFeo.25S in this study with those found in the body of existing literature presented in 


Table 5, we find consistent results [22]. 
3.1.3.4.2 Exchange Constants 


We have determined two important constants that shed light on the important role that the valence 
and conduction bands play in the spin-splitting exchange coupling observed in the electronic band 
structures of our compounds. These are the exchange constants for s-d and p-d, denoted Nog and Nog, 


respectively, which relate to CB and VB. Assuming the standard Kondo interaction, both constants were 


found [60]: 


AE, 

Noa = xess (eq.24) 
AEy 

Nog = 5s, (eq.25) 


In this context, AE, = E,\- E,| and AE, = E,'- E,', are the spin separation at the edges of the 
conduction and valence bands, respectively. The variable "x" denotes the percentage of Fe in the 


compound, while <S> signifies the average magnetization of the Fe atoms. 


Table 5 shows the calculated values for AE,, AE,, Nog, and Nog of the HMF compounds Srj.xFexS 
(x = 0.125, 0.25 and 0.50). Notably, the s-d exchange coupling between the CB and d-Fe states 
consistently exhibits an anti-ferromagnetic nature at all concentrations, as evidenced by the negative 
sign of the Nog constant. Conversely, the Nog constant shows a positive value for x = 0.125 and 0.50, 
but a negative value for x = 0.25. This means a mixture of ferromagnetic and antiferromagnetic 
interactions in the p-d exchange between the VB and d-Fe levels. The overall increase in Nogand 
fluctuations in Nog values with the increase in Fe content from 0.125 to 0.50 reinforce the targeted 


magnetic properties of these materials. 


Now 18 generally found to be positive in most studies on half-metallic ferromagnetic/semiconductor 
compounds, while Nog is invariably negative [44, 61]. In our study, this pattern is completely reversed. 
Changes in symmetry and quantum confinement effects may be responsible for the reversal of the sign 
of Nog [62]. In contrast, in the DMS Zn;.x Cr,Se, Mac and co-workers first predicted the positive sign 


of Nog and associated this significant shift with the d-state position of the valence band [63]. A 
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comparable pattern was also recently observed at low concentrations of the ternary compounds Cdo.937 


TMo.0625S (TM = Cr, V, Cu and Sc) [64]. 


Table 5 


The calculated values of total magnetic moment Mio (in UB) per Fe atom, local magnetic moments (in 


Us) of Fe (Mre), Sr (Ms;), and S (Ms), and magnetic moment in the interstitial sites M; (us), spin-splitting 


edges of conduction AE, (eV) and valence AE, (eV) bands, and exchange constants NO,, and NOg of 


Sris«FexS compounds (x = 0.125, 0.25, 0.50, and 0.75). 


Compound Met Mre Ms, Ms M; AE, AE, NO,  NOg 
Sros7sFeo.sS 4 3.690 0.005 0.031 0.119 -3.021 0.017 -12.08 0.068 
Sro.7s Feo.2s8 4; 3.680; 0.006; 0.052; 0.230; -2.639  -0.977_ -5.278 -1.954 

44 3.6512 — -0.0014_ —-0.044_—0.089 
Sro.so Feo.soS 4 3.643 0.010 0.125 0435 -2.768 0.884 -2.768 0.884 
Sro.2s Feo.7sS 3.973 3.636 0.004 0.158 0.750 : : : : 
‘Ref [22]. 


3.1.4 Conclusion 


In the first section of this chapter, we examined an in-depth study of the Fe-doped SrS 


semiconductor using rigorous first-principles calculations facilitated by the WIEN2K computational 


framework. This study involved a comprehensive investigation of the Sri.x.FexS compounds, where x 


takes values of 0, 0.125, 0.25, 0.50, and 0.75, with a focus on their behavior in the rock-salt phase. The 


key findings and conclusions from this endeavor are presented as follows: 


The results of the structural properties obtained with the PBE-GGA approximation showed that 
the ferromagnetic phase is the most stable in all Sri..Fe,S compounds and that the equilibrium 
lattice constant of the compounds decreases with increasing Fe concentration. The enthalpy of 
formation values confirmed the stability of the cubic structure, with Sro.soFeo.s0S being the most 
stable. 

Analysis of mechanical and thermal properties confirmed that all the doped-compounds meet 
the criteria for mechanical stability. They exhibit ductility and anisotropy and are characterized 
by ionic bonding. In addition, they exhibit remarkable resistance to heat treatment. Among 
them, Sro.so Feo.soS showed the highest stiffness, while Sro.g7sFeo.125S was found to be the most 


flexible. Furthermore, Sro.2sFeo.75S exhibited the highest degree of anisotropy. 
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3. The electronic properties calculated from the equilibrium lattice parameters confirm that all Sri. 
xFex,S compounds are HM ferromagnets and therefore applicable to spintronics, whereas the 


Sro.25Feo.75S compound shows metallic behavior. 


4. For all HMF compounds, we calculated a total magnetic moment of 4 pp and a low local 
magnetic moment at the non-magnetic sites of Sr and S due to hybridization between Fe-3d and 
S-3p states. The exchange interactions in the compounds yielded exchange constants Nog and 
Nog With opposite signs, confirming the existence of an opposite interaction between valence 


and conduction states. 
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3.2 Structural, Electronic, Magnetic, Optical, and Thermoelectric Properties of co-doped 
SrS: (Fe, p®) Alloys 


3.2.1 Introduction 


A key factor in the recent boom in scientific research has been the search for new materials with 
remarkable properties. In this regard, third-generation semiconductors — also referred to as — wide band 
gap (WBG) materials have emerged as a highly favorable platform for the development of high- 
performance optoelectronic and electronic devices [65]. SrS is a material that stands out among the 
others as an example with interesting properties. In addition to the features discussed previously in this 
Chapter, SrS is particularly known for its remarkable optical properties, which are covered in further 


detail in Chapter 1, Part 1.4. 


Despite extensive research on pure SrS, obstacles such as the large bandgap and indirect nature limit 
its application in optics [66]. Chemical doping is a promising solution to this problem, as demonstrated 
by our previous work, where doping transition metals, particularly Fe at 12.5%, successfully induced an 
indirect to direct band transition. This is a promising direction to improve the optical properties and 


performance of spintronics. 


The need for more data storage capacity makes it imperative to convert waste heat into usable 
electrical energy [67]. This plan is considered the best as it provides a reliable, cost-effective, clean and 
environmentally friendly source of renewable energy in the long term. High-performance spin- 
dependent thermoelectric materials lie at the intersection of thermoelectricity and spintronics, an area 


that requires further research and development. 


In particular, the alkali metal series has shown promise in doping the SrS site. These non-magnetic 
elements have been shown to increase ferromagnetism and luminescence efficiency [68-70] and are 
non-toxic. Simultaneous doping of these elements with transition metal ions in bulk SrS will be 


important as it can create new systems with unique properties. 


The second section of this chapter builds on our previous work and examines the structural, 
electronic, magnetic, optical, and transport properties of Li, Na, and K co-doped at a fixed concentration 
of 12.5% in SrS: Fe. Our knowledge of co-doped SrS systems is expected to grow as a result of this 
research, leading to the development and improvement of SrS-based materials for applications in 


thermoelectricity, magnetism, electronics, and optics. 
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3.2.2 Computational Details 


Since our work in this section of the chapter builds on the work in the previous section, we have 
chosen to continue with the same data parameter settings as for the mono-doped systems, including 
code, basis set, approximation, and function for exchange-correlation term within the INPUT 
parameters. We simply enlarged the supercell to accommodate 32 atoms in a (2 x 2 x 2) configuration 
of a bec structure in the FM case to produce the co-doped systems SrS: (Fe, p°= Li, Na and K). This 
supercell contains 16 Sr atoms and 16 S atoms. We exchanged two atoms at the Sr positions to achieve 
a total impurity ion concentration of 12.5%. One was fixed as Fe and is located at the (0, 0, 0) position 
(at the corner), and the other, for in ions, is located in the next position at (0.5, 0, 0), i.e. Sro.s7sFeo.0625 
p’o.062sS. For the AFM configuration, we used a larger P-type supercell with 64 atoms (2 x 2 x 2). To 
obtain trustworthy results in evaluating optical and thermoelectric properties, we also used a larger k- 
mesh and the TB-mBJ method [16]. We used Bardeen and Shockley's deformation potential theory [71] 
and Boltzmann's semi-classical theory-based BoltzTrap2 software [72] to perform a post-DFT approach 


to study the thermoelectric properties of the compounds. 
3.2.3 Results and Discussion 
3.2.3.1 Structural Properties 


The (2 x 2 x 2) bcc-lattice crystal of the (Fe, p®) co-doped SrS systems is modeled in Figures 11 (a) 
and (b). Table 6 shows the structural parameters details of the SrFeLiS, SrFeNaS, and SrFeKS systems, 
including lattice parameter (a), bulk modulus (B), total energy differences between FM and AFM 


configurations (AE=Earm -Erm), and Curie temperature (Tc). 


Figure 11. (a) Polyhedral view of the SrS: (Fe, p°) [p° = Li, Na, and K] crystal structures (b) The 


stable FM configuration. 
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Table 6 


The calculated values of crystal lattice parameter a (A), bulk modulus B (GPa), energy difference (AE), 
Curie temperature (T;), and formation energy (Ef) of (Fe, p°) co-doped SrS systems. The Fe mono- 


doped SrS data are included in the Table from the preceding section, to facilitate further analysis and 


comparison. 
Materials ay (A) Bo(GPa) = AEarmem (eV) ~—- Tc (KX) E; (eV) 
SrFeS 5.949 50.244 1.160 - -4.384 
SrFeLiS 5.983 48.322 0.432 1326 -5.099 
SrFeNaS 5.994 48.107 0.435 1337 -5.087 
SrFeKS 6.017 47.518 0.437 1341 -5.073 


The AFM and FM order’s energetic stability is indicated by the “-” and/or “+” signs of AE. According 
to the simulation results in Table 6, the positive AE values indicate that the co-doped SrS (Fe, p°) 
structures have the ability to suppress the AFM order and promote the FM order in the long term. This 
implies that the FM state is stabilized and gains an energetic advantage over the AFM state when Fe is 


combined with alkali metals (Li, Na and K). 


The Curie temperature (Tc), an important parameter for magnetic materials in spintronics, was 
calculated based on the energy difference JE. To estimate Tc, the MFA (Mean-Field Approximation) 


formula was used as follows: 


2 
TMFA — _ Inky OE (eq.26) 


Here kz is the Boltzmann constant and n is the number of substituted atoms. 


However, it is important to keep in mind that percolation behavior, which has a significant impact 
on the magnetic order and leads to an overestimation of T., is not taken into account by MFA [73]. To 
overcome this limitation, we used an empirical relationship [74] that uses the ratio T./T!*4= 0.794 to 
relate the exact value of the critical temperature (T,) to the expected mean field value abs 4) From 
Table 6, it is noted that all three systems have T. above room temperature (RT) and that there is a clear 
trend in T.. In particular, it increases as the atomic number (Z) increases from Li to Na to K, suggesting 
that the higher the Z, the more stable the system is. This highlights how the types of alkali metal dopants 
(Li, Na and K) can be changed to control ferromagnetism in these systems, opening possibilities for 
spintronic applications. It should be noted that this is the first time these systems have been predicted 


and further experimental validation is required. 
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Table 6 shows the optimized lattice parameters and bulk modulus of the co-doped systems in the 
stable FM configuration. Figure 12 shows the comparison of the total energies and the cell volume. The 
structural parameters of monodoped SrS: Fe** calculated in the previous section were also included to 
assess the feasibility of co-doping. Since the substituted Fe** (0.55 A) has a smaller ionic radius than the 
dopants Lit (0.76 A), Na* (1.02 A) and K* (1.38 A), the lattice constant increases gradually with 
increasing Z after the co-doping procedure. This indicates grid expansion. In contrast, the values of the 
bulk modulus decrease slightly, indicating that the co-doped systems have lower hardness. Importantly, 
Figure 12 shows that the addition of Li, Na, and K results in an increase in cell volume of 0.49%, 0.64%, 
and 0.98%, respectively, and a consistent energy reduction of about 3% to 8%. These data clearly 


indicate that the co-doping process can lead to significant stabilization of the system. 
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Figure 12. An illustration of the total energy and cell volume of the Fe-p° co-doped systems in 


comparison to the Fe mono-doped system. 


To evaluate the thermodynamic stability of the co-doped systems, we also calculated the formation 


energy. Our expression was as follows: 
Es = Ecodoped — AE” (Sr) — bE? (Fe) — cE?(p°) — dE?(S) (eq.27) 


The term Egodopea refers to the total energy of SrS co-doped with transition-alkali metals; the terms 
E®(Sr), E(S), E’(Fe), and E? (p°) denote the energies per atom of the host atoms (Sr and S), the 


fixed dopant (Fe), and the associated impurity (p° = Li, Na, and K); the numbers of the relevant atoms 


within the supercell are indicated by a, b, c, and d. 
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Table 6 shows that the estimated formation energies (EA of all compounds have negative values. 
This suggests that these substances may be synthesized under ambient conditions and are 
thermodynamically stable. Ey (re-Li) < Ey (Fe-Na) < Ey (Pre-K) < Ey (Fe-single) 18 the order in which the £ values 
are found. This implies that the most difficult system to dope is the singly Fe-doped SrS, which has the 
highest EZ On the other hand, the co-doped Fe-Li system has a comparatively simpler doping process 


as much less energy is required for its formation. 


In summary, we find that the introduction of monovalent alkali metals reduces the Ey of the mono- 
doped system. This implies that the alkali atoms help alleviate some of the difficulties in the singly Fe- 
doped SrS system, which in turn leads to a more stable configuration. Therefore, the overall stability of 


the system is significantly improved and the doping process is simplified by (Fe, p°) co-doping. 


3.2.3.2 Electronic Properties and Chemical Bonding Analysis 


3.2.3.2.1 Electronic Band Structure Analysis 


Figure 13 shows the electronic band structures of the co-doped SrS systems (Fe, Li), (Fe, Na) and 
(Fe, K) along the high symmetry directions (A, A, and Z) and high symmetry points ('-H-N-I°-P) of the 
first Brillouin zone (BZ). 


All SrS co-doped systems are half-semiconductors because the Er is clearly empty since neither the 
up- nor the down-spin channels crosses it. Both the VBM and CBM are positioned at the J" point for up- 
spin, resulting in a wide direct bandgap of 3.241 eV, 3.246 eV, and 3.182 eV for co-doping of Li, Na, 
and K, respectively. On the other hand, the narrower and indirect bandgap is indicated by the CBM 
being at the H point and the VBM being at the J’ point for the down-spin channel. For co-doping with 
Li, Na, and K, the bandgap size upon spin-down was found to decrease with Z, with values of 0.353 eV, 
0.346 eV, and 0.300 eV, respectively. It is important to note that the VBM is located close to the Er in 
the band structures of all three compounds, suggesting that they are p-type semiconductors. The 
comparatively smaller bandgap of the spin-down channel (<1 eV) facilitates thermal excitation of 
electrons, which improves electrical conductivity. Due to this property, these p-type HSC alloys have 


great potential for use in thermoelectric (TE) systems. 


Table 7 summarizes the results for the up-spin and dn-spin channels. We find that the introduction 
of hole-type charge carriers through alkali doping not only maintains the semiconducting nature of the 
up-spin but also reduces its bandgap values compared to Fe-doped SrS, which exhibits HMF properties. 


It also causes a down-spin shift from the Er, leading to conventional semiconductor behavior. 
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Table 7 


The calculated bandgap values of co-doped SrS: (Fe, p®) systems. The data of Fe mono-doped SrS are 


included in the table from the previous section to facilitate further analysis and comparison. 


Materials Gap Up-spin channel Down-spin channel 
nature Direction Value (eV) Direction Value (eV) 
SrFeS HMF (l-T] 3.401 = = 
SrFeLiS HSC (T-T] 3.241 [I'- H] 0.353 
SrFeNaS HSC (l-T] 3.246 [I- H] 0.346 
SrFeKS HSC [T-T] 3.182 [I'- H] 0.300 


Energy (eV) 


aw 


Energy (eV) 


Energy (eV) 
Energy (eV) 
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Figure 13. (a-c) Spin-resolved electronic band structures of (Fe, Li), (Fe, Na), and (Fe, K) co-doped 
SrS with mBJ-PBE method. 


3.2.3.2.2 Total and Partial Densities of States Analysis 


The spin-resolved total density of states (TDOS) and orbital-decomposed density of states (PDOS) 


of Sr, S, Fe, and p” atoms are shown in Figure 14. 


DOS (state/eV) 
DOS (state/eV) 
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Figure 14. (a-c) Spin-resolved total and partial densities of states of of (Fe, Li), (Fe, Na), and (Fe, K) 
co-doped SrS with the mBJ-PBE method. 


The asymmetry of the TDOS in both spin channels can be seen in Figures 14(a—c) and confirms the 
ferromagnetic nature of the materials. It can also be seen that neither spin channel exhibits an impurity 
state at the Er level, highlighting the magnetic semiconductor nature of them. Fe:3d and S:3p orbitals 
make up the majority of the VB, while the S:3p and Fe:3d orbitals dominate the VBM, with the Sr:5s 
and p®:s (Li:2s, Na:3s and K:4s) orbitals contribute very little. Fe:3d orbitals make up the majority of 
the CB, and Fe:3d determines the CBM, witha small contribution from Sr: 5s and p®:s orbitals. It appears 
that alkali substitution contributes only slightly to the TDOS near the Er. Ultimately, this is different 
from the TDOS detected by the Fe mono-doped SrS, where the compound exhibited the HMF property. 
We hypothesize that coupled Fe-p° doping is preferable to Fe mono-doping to achieve optimal magnetic 


and optical properties due to the presence of a localized energy level at the top of the CB near the Ep. 


A careful examination of the PDOS shows that the Fe:3d values have shifted significantly compared 
to the parent compound SrS:Fe. Atomic Fe:3d levels split into eg and t2g manifolds with lower and 
higher energies, respectively, due to the tetrahedral crystal field. While the e, orbital splits into the 
singlet d,.2 —y2 and d,2 levels, the tz, orbital splits further into the double degenerate dx, + dy, levels and 
the singlet dxy level. Hund's rule states that since the 3d state of the Fe dopant has six electrons, the five 


up-spin electrons fill the five sub-states that lie below the Er level in a potential range of about ~-0.4 to 
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-4.5 eV. Since the p®-s states are one electron short, the addition of alkali dopants creates a hole in the 
system. As a result, the atoms of Li, Na and K act as acceptors. Thus, the p-type semiconductor behavior 
in the three compounds can be explained by the sub-states in the spin-down channels remaining 
unoccupied and located above the Er level in a potential range of approximately ~ — 0.3 eV—1.4 eV. 
The 3p orbitals of neighboring S atoms and the 3d orbitals of Fe clearly overlap close to the Ep, 
suggesting strong pd-hybridization between them. This pd-hybridization is crucial for maintaining 
ferromagnetic order and generating magnetism in the co-doped systems in addition to band edge 
splitting. Another factor is the coupling chain between the states Fe:3d,2¢/g—Sr:5s—p":3s, which exists 


in the energy range from about ~-2 eV to Er. 


Along with the previously defined direct exchange splitting (A,(d)), we also calculated the crystal 
field energy (AE;;-y,), which is defined as the energy arising from the interaction of surrounding ligands 
with metal ions (E,2, — Egg). When the A,(d) energy exceeds the AE,,-y,, the energy obtained from the 
parallel spin alignment through direct exchange is larger than the energy difference between the e, and 
tz, orbitals, this can cause the material to exhibit more pronounced ferromagnetic behavior. In our case, 
SrFeLiS, SrFeNaS and SrFeKS have calculated A,(d) values of 6.06 eV, 6.14 eV and 5.465 eV, 
respectively. These values exceed their corresponding AE;;y; values of 0.672 eV, 0.713 eV and 0.436 
eV. Therefore, they strongly support the parallel alignment of the spins and promote the stabilization of 


ferromagnetism in these systems. 


A crucial parameter for understanding the electronic structure and optical properties of (Fe, p°) co- 
doped SrS is the d-d transition, which is defined as the distance between up-spin occupied Fe:3d states 
and unoccupied Fe:3d states in down-spin. We gain important insights into the energy absorbed during 
the d-d transition by quantifying the energy separation between the up-spin and down-spin states 
involved in this transition. Our calculations show that the co-doped SrS systems (Fe, Li), (Fe, Na), and 
(Fe, K) have d-d transition energies of 0.300 eV, 0.346 eV, and 0.353 eV, respectively. These values 
show that the co-doped SrS systems absorb 0.300 eV, 0.346 eV, and 0.353 eV during the d-d optical 


transition in the FM configuration, respectively. 
3.2.3.2.3 Electron Density Difference and Chemical Bonding Analysis 


The analysis of the electron density difference of the co-doped systems provides information about 
the interaction between atoms and their bonding properties. In an ionic bond, there is a significant 
difference in electronegativity between the bonded atoms, resulting in electron transfer and the 
formation of charged ions. This type of bond is predominantly ionic in nature. Conversely, a covalent 
bond occurs when atoms with similar electronegativity share electrons, resulting in a more balanced 


electron distribution and covalent character. 
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The 2D bond charge density distributions of the SrS: (Fe, Li), SrS: (Fe, Na), and SrS: (Fe, K) systems 
were plotted in the (110) plane in Figure 15. The charge density difference is represented by the color 


scale in the figures, with red indicating electron depletion and blue indicating electron accumulation. 
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Figure 15. (a-c) The bonding contour profiles of electron density difference for (Fe, Li), (Fe, Na), and 


(Fe, K) co-doped SrS Systems, respectively. 
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A small amount of electrons, about 0.016 eV/A3 accumulates between Sr and S, Fe and S, and p” 
and S, as shown in this figure. This indicates the strong ionic character of the Sr-S, Fe-S and p°-S bonds. 
Interestingly, there is a larger accumulation of electrons around the Fe atom compared to the Sr atom, 
suggesting that Fe has a greater ability to donate or accept electrons, which affects its bonding behavior. 
Furthermore, the analysis showed that more charges were localized between p° and S, reducing electron 
exchange between them. As a result, there is a greater concentration of free electrons generated by p” 
atoms, leading to an increase in electron accumulation near p°. For co-doping with Li, Na, and K, the 
electron accumulation values around p” are about 0.05, 0.196, and 0.3 eV/A3, respectively. This suggests 
that the ionic character of the p°-S bond in the (Fe, p®) co-doped SrS system slightly improves with 
increasing Z of alkali metals. The doping of p° atoms thus improves ionic bonding of SrS: (Fe, p°) co- 


doped systems. 
3.2.3.3 Magnetic Properties 
3.2.3.3.1 Magnetic Moment 


Table 8 summarizes the results of the total magnetic moment, the magnetic moments of Sr, S, Fe 
and p°, and the magnetic moment of the interstitial region calculated by the PBE + mBJ method. The 
calculated Mror in all systems converges to an integral value of 5.00 Ls, as shown in this table. The 
mono-doped SrS:Fe indicates a permanent magnetic moment of 4.00 [1p in free space, which is different 
from this value. As mentioned in the DOS section, the total value of the magnetic moment increases to 
5.00 ps due to the charge exchange between the Fe and p’ ions in the co-doped systems. The 3d orbitals 
of the transition metal ion (Fe) and the interstitial region are the main sources of Mror. The nearest 


neighboring S and Sr atoms as well as the alkali metals contribute very little to the total value. 


The sp-d hybridization between the Fe:3d states and the S:3p, Sr:5s, and p®:ns states (where n is the 
number of Li, Na, and K electrons) in these compounds is responsible for the higher values of magnetic 
moments in the interstitial region. Due to this hybridization, the magnetic moments of the interstitial 
region increase and the local magnetic moments associated with the Fe atoms decrease (to about 3.89 
up in all compounds). This phenomenon illustrates how hybridization has a significant impact on how 


these compounds behave magnetically and distribute electrons. 


The studied co-doped systems exhibit elevated and integer magnetic moment values, suggesting 


possible applications in spintronic devices. 
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3.2.3.3.2 Exchange Constants 


The s-d and p-d exchange constants are determined using the corresponding equations (Eq. 24) (Eq. 
25) described in the first section of this chapter. The calculated AE,, AEy, Nog, and Nog are summarized 
in Table 8. The negative values of s-d coupling suggest an attractive interaction between magnetic 
moments on dopant atoms and band holes, implying AFM coupling in the alkali-doped SrS:Fe systems. 
Conversely, the positive values of Nog indicate that the p-d coupling is repulsive, indicating FM 
behavior. Compared to the results of Fe single doping, the exchange constants for p’co-doping show a 
slight increase. This slight increase can be attributed to the effects of additional electrons, altered 
hybridization, and variations in the crystal lattice and local environment. These factors collectively 


influence the magnetic coupling and contribute to the observed changes in exchange constants. 


Table 8 


The calculated values of total magnetic moment (M™), interstitial magnetic moment (M'), magnetic 
moment on each Sr, S, Fe, and p° atoms, band edges splitting AE, and AE,, and exchange constants 
Noqwand Nog of co-doped SrS systems. The data of Fe mono-doped SrS are included in the Table from 


the previous section for further analysis and comparison. 


Materials Mt (up) Mi (us) M¥e(us) MS*(up) MS(us) Mus) AE-(eV) AE, (eV) NO, NO, 


SrFeS 4 0.119 3.690 0.005 0.031 - -3.021 0.017 -12.08 0.068 
SrFeLiS Bs) 0.503 3.893 0.008 0.110 0.002 -2.841 0.053 -8.782 0.170 
SrFeNaS 8) 0.270 3.895 0.007 0.111 0.112 -2.840 0.052 -9.088 0.171 
SrFeKS re) 0.268 3.892 0.007 0.112 0.005 -2.830 0.048 -9.05 0.154 


3.2.3.4 Optical Properties 


The optical behavior of a material in a medium is elucidated by examining its optical properties. 
Highly spin-polarized semiconductors (HSCs), which span a wavelength spectrum from ultraviolet to 
visible, can be treated as a continuous medium. Under optical excitation, these materials demonstrate 
the ability to generate fully spin-polarized electrons and/or holes. The complex dielectric function €(@) 
is a key optical property that bridges the gap between the physical process of electron transfer and the 
energy band structure. It reflects the linear macroscopic response of the material to electromagnetic 


radiation. The complex dielectric function can be expressed as: 
E(w) = €1(w) + i€2(w) (eq.27) 


Where €2(w) signifies the imaginary part of the complex dielectric function and delineates optical 


absorption. Conversely, the real part €;(@) denotes dispersion and polarized radiation. 
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The €2(w) includes two different contributions: interband transitions (between the bands) and 


intraband transitions (within the bands). It can be calculated based on the band structure as follows [77]: 


An” e2 


E2(@) = (“F-5) Dig S(UMIP)? Fi — Fp)S(Ep — Ei — w) dP k (eq.28) 


m2 w2 


Where f; and f;, are the Fermi functions for initial |i) and final |f) states with energies E; and Ey and 


M is the momentum operator. 


The €;(q@) can be derived from the imaginary part through the Kramers-Kronig relation [75], as 


shown below: 


€\(w) =1+ =p ie wre) aay! (eq.29) 


0 w-qw2 
Where P is the Cauchy principal value. 


As @ > 0, the static dielectric constant (€;(0)) in Eq.29 simplifies to: 


€1(0) = e(w) =1+=P fr 2@ do’ (eq.30) 


The integral part dee 200) doo! is proportional to a value of + Wp. The @, is defined as the plasmon 


frequency [75] and it is given as follows: 


2 Nye? 


w (eq.31) 


P MEQ 
Where N,,, is the valance electron density and m@ is the effective mass of the electrons. 


From the components €;(@) and £€2(@) of the complex dielectric function, further parameters 


such as the refractive index n(w) and the extinction coefficient k(w) could be calculated [75]: 


L e3(w) +8 (w)+e4(w)]'/2 

n(@) = = (eq.32) 
[ | e2 (co) +22 (wo) £4 (co)] 2 

k(w) = >_> (eq.33) 


v2 


We can also calculate the absorption coefficient a@(w) and the optical conductivity o(@) using the 


following formulas [75]: 


a(w) = V2 (w)[ ES (w) + €3(w) — €1(w)] */2 (eq.34) 
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o(w) = A e,(w) (eq.35) 


Figure 16 shows the real and imaginary parts of the complex dielectric function, the absorption 
coefficient and the optical conductivity, over a photon energy range of 0—13.6 eV. In addition, for 
comparison purposes, we calculated the optical properties of pure SrS and Fe-doped SrS using the data 


from the first part of this chapter. 
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Figure 16. (a-d) Real part of the complex dielectric function €;(@), Imaginary part of the complex 
dielectric function €2(@), Absorption coefficient a(@), and Optical conductivity o(@) spectra, 


respectively, of pristine SrS, Fe-single doped and (Fe, p°) co-doped systems. 
a. Real and Imaginary Parts of Complex Dielectric Function 


The discussion begins with an examination of the real part of the complex dielectric function. In 
Figure 16(a) we plotted the profiles of the ¢:(@) function for the binary, ternary, and quaternary systems. 


At lower energies, a notable difference is apparent between the co-doped, singly doped, and pristine SrS 
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systems. The static dielectric, denoted by €1(0), refers to the real part value at zero energy and is defined 
by (Eq. 30). In the case of the pristine SrS, ¢:(0) starts positive and reaches a value of 3.769. After 
reaching this value, €:(@) increases with energy until it reaches its maximum of 6.151 at 4.186 eV in the 


visible region. 


When Fe is introduced into pure SrS, the compound shows the HME character, which causes €)(0) 
to shift towards a significantly large negative value, indicating two crucial aspects: First, this negative 
value means metallic behavior in this energy sector. Secondly, it results in energy loss and limited light 
transmission through the medium, resulting in significant reflections. However, once it reaches its 
minimum value, the magnitude of €:(@) begins to increase and finally becomes positive at 0.416 eV, 


marking the first root. 


When co-doping SrS:Fe with alkali metals, sharp peaks appear on the low energy side, showing 
high values of €:(0) with 6.194, 5.954 and 5.822 for Li-, Na- and K-co-doped SrS:Fe, respectively. These 
values are significantly increased compared to pure SrS. These distinct peaks indicate the 
semiconducting behavior that arises upon the introduction of alkali metal ions. The significantly higher 
value of €:(0) in SrFeLiS underlines its superior polarizability compared to other alloys. Beyond 0.416 
eV, the €:(@) profiles of the three co-doped systems follow a similar trend as the pure SrS and Fe single 
doping. The maximum dielectric peak for Fe mono-doping occurs at a photon energy of 4.3 eV, which 


corresponds to a dielectric value of 6.204. 


Beyond this energy, there is no apparent disparity in the €1(@) profiles of all studied systems. This 
suggests that the influence of alkali metal co-doping may also extend to lower energies. Subsequently, 
the value of ¢:(@) decreases sharply and finally reaches negative values after 8.065 eV. These negative 
values imply high reflectivity and plasmonic excitation in the energy range of 8.065—13.6 eV, indicating 
metallic behavior of the materials. Furthermore, the €:(@) values show a trend towards greater stability 
with increasing photon energy. This suggests that the materials’ interaction with incoming high-energy 


photons becomes more restricted or dampened. 


Figure 16(b) provides a visual representation of the spectra of the €2() function. Such spectral data 
provide crucial insights into how materials react to incoming electromagnetic radiation. In addition, the 
phenomenon of interband shift is clearly illustrated, which arises from the interaction between electrons 


in the occupied valence band and the radiation built into the material below the Er. 


In contrast to the behavior of €1(@), €2(@) in the Fe mono-doping system shows high positive values, 
followed by a sharp decay, with multiple transition peaks appearing at different rates in the low energy 
region. These interband transitions can be attributed to the influence of impurity bands introduced by 


the electrons of the Fe atom. The decrease in €2(@) indicates that the interaction between photons and 
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the material surface is minimal [76]. Prominent peaks in the €2(@) spectrum of Fe-doped SrS indicate 
electronic transitions between occupied and unoccupied states. These transitions can be identified 
directly by the PDOS of each connection. The 2.0 eV peak arises from the transition between 3p-S in 
the VB and 3d-t2, Fe in the CB. At 2.7 eV, an electronic transition occurs between the 5s-Sr state and 
the 3d-t2, Fe state. The significant 5.4 eV peak is attributed to the electronic transition between the 3d- 
t2z orbital in the VB and the 3d-e, orbital in the CB (see Figure 9 in this chapter). The variation in the 
size of these dielectric peaks can be explained by the difference in the amount of photoelectron energy 


required for electronic transitions induced by the different dopants. 


When examining the €2(@) curve of pristine SrS, a prominent peak is visible in the near UV region 
and maximum absorption occurs in this region. The absorption edge for this binary system begins at an 
energy of 3.435 eV and is due to electron transitions between the S-3p states in the VB and the Sr-5s 


states in the CB. 


For the systems SrFeLiS, SrFeNaS and SrFeKS they show clear absorption peaks in the infrared 
range of 0.3-1.6 eV. This means enhanced infrared absorption in the SrS and SrS:Fe systems through 
co-doping with one of the three alkali metal ions. It is noteworthy that co-doping with Li has the highest 
intensity. The absorption edges in these ordered alloys SrFeLiS, SrFeNaS, and SrFeKS manifest as 
energy peaks at 0.353 eV, 0.346 eV, and 0.300 eV, respectively. These peaks arise from electron 
transitions between the S-3p states in the VB and the Fe-3de, states in the CB. Remarkably, the positions 
of these energy peaks closely match the forbidden gaps of the respective alloys. This close agreement 
suggests that the electron transitions responsible for the absorption edges occur within the energy range 


defined by the forbidden gap of each material. 


sp-d hybridization creates direct transitions from the VB to the CB, leading to the formation of 
shallow levels near the Er. This reduces the bandgaps of the materials, which leads to a red shift. In the 
energy range of 1.6—2.7 eV, the co-doped systems exhibit higher peak intensities compared to the pure 
and singly doped systems, which means that the co-doping enhances the absorption of the system in the 
visible region. Significant peak intensities of €2(@) are observed in the energy range of 3.2—11.3 eV, 
indicating significant absorption from the near visible to ultraviolet region. This property makes the 


studied materials promising for applications in energy production systems. 


b. Absorption Coefficient 


In addition to €;(@) and €2(@), we calculated the absorption coefficient o(@), which is a crucial 
metric that quantifies the material's ability to absorb and convert photons of specific energy into a usable 


form. This coefficient provides valuable insight into the attenuation of incident light intensity as it passes 
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through a material and illustrates the process of light decay in the absorbing medium. The a(@) curves 


of pristine SrS, Fe-doped, and (Fe, p°) co-doped SrS are shown in Figure 16(c). 


The HMF property is evident in the absorption coefficient spectra. In the singly doped system there 
is a clear infrared absorption in the lower energy range from 0 to approximately 0.4 eV. However, a 
notable absence of infrared absorption is observed in the pure and co-doped systems. This is because 
the photon energy falls into the forbidden bandgap, which means transparency for infrared light. As the 
photon energy increases, the absorption coefficient increases rapidly, showing typical semiconductor 
behavior. Of all the compounds examined, the ultraviolet range has the highest absorption. The 
absorption spectra of these materials show prominent peaks in the UV range, which mainly occur above 
9 eV photon energy. These peaks show specific values, namely 162, 163, 162, 159, and 157 [10*Cm—1], 
for the pristine Fe-mono-doped and Li, Na, and K co-doped SrS alloys, respectively. This indicates a 
significant level of absorption and minimal electron losses within this energy range for all alloys. In the 
visible range, the absorption values for all compounds are comparatively lower than in the UV range 
and almost zero for SrS. The substitution of p° also shifts the position of the absorption peaks of singly 
doped SrS in the infrared and visible region and results in new peaks. Consequently, a redshift occurs, 


which is consistent with the previous interpretation of €2(@). 


Based on previously reported data, the absorption curves can be divided into two distinct regions 
showing improvement compared to pure SrS. SrS: (Fe, p°) compounds exhibit significant absorption 
compared to SrS: Fe in the first region, covering the energy range from 0 to 2.7 eV and covering the 
near-infrared to visible region. In the second region, which extends from 2.7 to 13.6 eV and corresponds 
to the near ultraviolet region, SrS:Fe absorbs more energy than SrS:(Fe, p°). These groundbreaking 
findings clearly demonstrate the potential of SrS: (Fe, p°) materials as highly effective hole transport 
materials (HTM) for use in solar cells. They are particularly characterized by their exceptional property 
of minimal absorption of incident sunlight, making them an ideal choice for improving solar cells 


efficiency. It should also be mentioned that SrS:Fe is ideal for use as a photodetector. 


c. Optical Conductivity 


The optical conductivity o(@) is related to the conduction of electrons, which are formed when a 
photon of a certain frequency interacts with a medium. Figure 16(d) shows the optical conductivities of 
pure, Fe-doped and (Fe, p°) co-doped SrS. The optical conductivity spectra largely agree with the shape 
of the imaginary part of the dielectric function. This means that by analyzing these spectra, which come 
from the imaginary part of the dielectric function, we can identify specific transitions that are responsible 


for the peaks. This is done using band-to-band decomposition as explained earlier. 
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6(@) of the co-doped compounds is quite intriguing and starts around the fundamental bandgap. 
Notably, the optical conductivity for all compounds, including the pure and singly doped ones, shows 
significant peaks in the photon energy range of 4-11 eV. This indicates an increased photon absorption 
rate in this particular energy range. Remarkably, all compounds in this region exhibit similar behavior, 
highlighting their comparable optical properties. The remarkably high value of optical conductivity 
means that these materials respond strongly to light, demonstrating their efficiency in conducting 
electrons when illuminated. Furthermore, the optical conductivity patterns consistently mirror those 
observed at o(@), highlighting that the absorption of photons results in an abundance of charges available 


for conduction. This highlights the crucial role of photon absorption in enabling charge conduction. 


It is important to highlight that our optical findings for the binary system agree well with previous 
theoretical evaluations [20, 77, 78]. However, it is worth noting that both experimental and theoretical 
data are lacking for other ternary and quaternary systems, highlighting the need for experimental 


validation. 
3.2.3.5 Thermoelectric Properties 


As stated in Section 1.3.3 of Chapter 1 and 2.3.3 of Chapter 2, thermoelectricity is the process of 
converting heat into electricity through the Seebeck effect. In this scenario, it is the heat flow that induces 
movement of charge carriers, resulting in the generation of an electric current. The thermoelectric 
performance of materials can be evaluated using the dimensionless figure of merit, which is described 


as follows: 
ZT = S’oT/(ke + Kl) (eq.36) 


In this case, § stands for the Seebeck coefficient, o for the electrical conductivity, PF = So for the 
power factor, Ke+ ki for the total thermal conductivity (which adds the contributions of the lattice and 


electronic components, respectively) and T for the temperature in Kelvin. 


To achieve the best thermoelectric performance, low thermal conductivity, significant Seebeck 
coefficient and high electrical conductivity are required. Nevertheless, the electronic component of 
thermal conductivity is directly related to electrical conductivity according to the Wiedemann-Franz law 
[79], which establishes a connection between ke and electrical conductivity via ke = LoT (where L is 
the Lorentz number). Because of this interdependence, it is difficult to improve ZT by changing one 


parameter without also changing the others. 


Selecting and designing thermoelectric materials, therefore, requires considering various other 
crucial factors. This includes aiming for a low band effective mass, which leads to a high Seebeck 


coefficient, without significantly impeding electrical conductivity. In BoltzTrap2's constant relaxation 
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time approximation (RTA), both electrical and electronic thermal conductivities, unlike the Seebeck 
coefficient, are dependent on relaxation time (7). Therefore, determining t is necessary to determine 
their actual values. The electronic structure derived from DFT calculations serves as an important input 
for the estimation of this parameter. Based on the identified band structures, we determined the effective 
masses (m*) of the charge carriers at the CBM and VBM using the parabolic band approximation 
framework [80]. Additionally, we calculated the uniaxial deformation potentials of electrons and holes 


by evaluating the CBM and VBM changes. 


We also performed estimates of uniaxial deformation potentials for electrons and holes obtained by 
calculating changes in the CBM and VBM. The deformation potential (DP), defined as the energy 
change of a system due to a deformation of its lattice structure along the transport direction, is described 


as follows: 


Cr) e e 
E,= te (eq.37) 


Where OE edge represents the alteration in CBM and VBM for electrons and holes, respectively, and < 
0 


denotes the uniaxial strain, ranging from -4% to 4%. The slope of the curves represents Ea. 


Based on their effective masses and strain potentials, we calculated the relaxation time, for the 
pristine, doped and co-doped SrS systems using the deformation potential theory proposed by Bardeen 


and Shockley [71]. This theory provides a formula that relates t to temperature as follows: 


_ av 2mh* Cy 
C= Sa sanz (mk, T)/2E (eq.38) 


h is the reduced Planck’s constant, m* is the band effective mass, Ea is the DP, and C, is the elastic 


constant computed via a quadratic polynomial fit of total energy fluctuation under strain as: 


2 
Ce (eq.39) 


~ 4 Aaya 
Vo Oo 
Herein, Vg represents the equilibrium volume. 


For the pure, doped, and co-doped SrS systems, the fitted curves showing edges and total energies 
versus strains are shown in Figures 17—20. Table 9 additionally shows the calculated values of m*, Ea 


and C,. 
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Figure 18. Band energy edge (Eeage) as a function of uniaxial strain (“) for: (a) Fe mono-doped and 
0 


(b-d) (Fe, Li), (Fe, Na), and (Fe, K) co-doped systems, respectively. 
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Figure 19. Total energy as a function of uniaxial strain () for the pristine SrS. 
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Figure 20. Total energy as a function of uniaxial strain (“) for: (a) Fe mono-doped and (b-d) (Fe, 
0 


Li), (Fe, Na), and (Fe, K) co-doped systems, respectively. 


As shown in Table 9, the Cy values gradually decrease as dopants are introduced into SrS and as 
the system goes from mono-doping to co-doping. This is consistent with an increase in atomic number 
within the co-doped systems. This decrease in the elastic constant can be attributed to the introduction 


of additional dopant atoms, which disrupt the regular lattice structure and lead to lattice defects. 


On the other hand, the Eq has a stronger effect on electrons than on holes. These different effects 


highlight the different sensitivity of these charge carriers to external perturbations. 


In addition, the m* of the majority-electrons is smaller than that of the majority-holes. This contrast 


in m* highlights the different mobilities and transport behavior of the majority electrons and holes. 


Comparing our results for the pristine material with those of Hou et al. using HSE [80], we observe 
good agreement for all parameters. The slight differences could possibly be due to differences in the 


computational code or approximation methods used (VASP code). 
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Table 9 


The calculated elastic constant Cg, the DP constant Ea, and the effective mass m* of the pristine, doped 


and co-doped SrS systems. 


Materials Carrier type Cy (eV/A3) Ea(eV) m* (me) 
SrS p-type 1.335 6.25 0.461 
1.3473 -17.203 0.590! 

n-type 1.335 -18.55 0.170 

1.3473 -11.44) 1.518) 


Up-Spin Dn-Spin Up-Spin Dn-Spin 


SrFeS p-type 0.800 4.951 — -1.150 0.538 1.363 
n-type 0.800 -12.25 -9.951 0.288 0.899 
SrFeLiS p-type 0.780 -3.452 -0.952 0.560 0.382 
n-type 0.780 -11.05 -3.703 0.254 0.841 
SrFeNaS p-typ 0.696 -1.080  -0.544 0.583 0.393 
n-type 0.696 -9.853 3.501 0.244 0.910 
SrFeKS p-type 0.644 0.734. --0.230 0.600 0.417 
n-type 0.644 -5.000 -3.100 0.227. 0.795 
i Ref [82]. 


Using the obtained parameters from Table 9 and substituting them into Eq. 38, we can determine T. 
The t of holes (p-type) and electrons (n-type) for pure SrS as well as for majority- and minority- spins 
of doped and co-doped systems as a function of temperature in a range of 100-1200 K is shown in 


Figures 21 and 22, respectively. 
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Figure 21. Measured relaxation time (t) of holes and electrons as temperature dependence for SrS. 
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Figure 22. Measured spin-resolved relaxation time (t) of holes and electrons as temperature 


dependence for (a) Fe-single doped and (b-d) (Fe, p°) co-doped systems, respectively. 


From the figures, it can be seen that all compounds have similar variations in t. The intensity of t 
initially decreases before gradually varying over a high temperature range. At low temperatures, the co- 
doped system SrFeKS has the most significant value, followed by SrFeNaS, SrFeLiS, SrFeS and then 
SrS. This trend can largely be attributed to the low deformation potential values of the SrFeKS material 
compared to other materials (since Tt is inversely proportional to the square of Eu, the impact of this 
contribution will be particularly significant). The relaxation time of HSC compounds is ten times longer 
than that of HMF compound. Consequently, we expect high electrical conductivity and consequently 


high ZT in all p-type co-doped systems. 


Our calculated t values for pristine SrS appear to be higher than those of Hou et al. [80], who 
obtained values in the order of 104s. This is also due to the difference in the computational code or 


approximation methods used 
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After determining t, we calculated the spin-dependent transport properties, including o, S, «, and 
ZT, for both holes and electrons in all systems, independent of t. The results are shown in Figure 23 (a- 


e). 
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Figure.23 Temperature dependent thermoelectric properties (a) electrical conductivity (0), (b) 
Seebeck coefficient (S), (c) electronic thermal conductivity (Kk), (d) lattice thermal conductivity (1), 
and (e) figure of merit (ZT) for both p-type and n-type of pristine, Fe-single doped and (Fe, p”) co- 
doped SrS systems. 


a. Electrical Conductivity 


Electrical conductivity (o) refers to the movement of free electrons within materials. It generally 
depends on the mobility of the charge carriers, which is influenced by impurities and defects. Electrons 
move from the hot to the cold region in a material, and as a result of this phenomenon, electricity is 
generated. For an effective thermoelectric device, materials must have high electrical conductivity. 
Materials are categorized as either conductors or semiconductors according to energy band theory. 
Conductors have freely moving electrons for conduction, while semiconductors require the supply of 


external energy to enable charge movement. Our analysis uses the two-current model to study the 
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temperature-dependent changes in spin-based transport properties. The plots of o for p-type and n-type 


are shown in Figure 23(a) versus temperature up to 1200 K. 


In the two-current model, electrical conduction is determined by taking into account the currents in 
both spin channels [81]. This means that the total electrical conductivity is the sum of the respective 


values for each spin channel, expressed by the equation: 6 = of + o| [82]. 


Examining the curves, it is clear that the electronic conductivity of pure material is almost 
independent of temperature for both p-type and n-type carriers. This observation is consistent with the 
findings of Hou et al. [80] and Rajput et al. [83]. Conversely, the curves for the co-doped systems all 
show a gradual increase with increasing temperature, which is characteristic of semiconductors. On the 
other hand, the conductivity of mono-Fe-doped SrS gradually decreases with increasing temperature, a 
behavior typical for metallic materials. Notably, the conductivity of p-type compounds exceeds that of 
the corresponding n-type compounds, which may be due to the longer relaxation time in p-type systems 
compared to n-type systems. For p-type carriers, the maximum conductivity reaches 103.1x10’ S.m'! 
for SrS: (Fe, K), followed by 7.4x10’ S.m'! for SrS: (Fe, Li), then 0.7x10’ S.m™! for SrS: (Fe, Na), 
1.6x10’ S.m' for SrS: Fe and a minimum of 0.0125x10’ S.m"! for pristine SrS, all at 1200 K. 


b. Seebeck Coefficient 


The Seebeck coefficient (S), commonly known as thermo-power, is quantified as the voltage 
produced in an open circuit between two points where there is a uniform temperature difference of 1 K. 
The magnitude of this coefficient varies with the temperature level at which the thermal gradient is 
introduced. An essential Seebeck coefficient is essential for optimizing thermoelectric devices. In the 
two-current model, the total Seebeck coefficient can be understood as the average of the up-spin and 
down-spin coefficients, weighted by their respective electrical conductivities as S = of St + o| S|/of 


+6] [82]. The variation of the Seebeck coefficient with temperature is shown in Figure 23(b). 


The sign of the Seebeck coefficient basically determines the type of majority charge carriers in the 
system. Figure 23(b) shows that both types of carriers exist in our considered series of materials: (i) 
materials with a positive Seebeck coefficient region have holes as the majority charge carriers, and (ii) 
materials in the negative Seebeck coefficient region have electrons as the majority of charge carriers. 
The Seebeck coefficients of p-type and n-type SrS, SrS: Fe and SrS: (Fe, p°) alloys show similar trends. 
In the range of 200-500 K, the absolute values of S for both carrier types in pure SrS initially decrease 
with increasing temperature, reaching a minimum of 245 pV K™ at 500 K and then decreasing with 
temperature until they reach a maximum of 249 reach uV K' at 1200 K, which is comparable to previous 
results [80, 83]. For the HMF compound, S increases gradually, reaches a maximum of 24.29 pV. K"! 


at 500 K, and then slowly decreases at temperatures between 500 and 1200 K. Based on these results, 
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we can conclude that doping iron in SrS is not preferable as it reduces its Seebeck coefficient, showing 
that pure semiconductors have better voltage generation compared to HMF compounds. For Li- and Na- 
co-doped systems, S increases with increasing temperature, indicating a high Seebeck coefficient at high 
temperatures. Conversely, co-doping with K initially decreases S in the range of 200-500 K, followed 
by a linear decrease to match the trend of pristine and other co-dopants. The simultaneous doping of Li, 
Na and K in SrS:Fe significantly increases the absolute value of S, approaching 151.15, 160.29 and 
150.24 nVK"! at 1200 K, respectively. The variation of S values can be attributed to the differences in 


the bandgap values. 
c. Thermal Conductivity 


Optimization of thermoelectric materials also depends on a crucial factor known as thermal 
conductivity (x). This factor refers to the transfer of energy in the form of heat due to temperature 
differences within a material. The total thermal conductivity (k= «ke + «l) includes both the electronic 
(ke) and the lattice contribution (x1). As already described, «e is closely linked to electronic conductivity 
via the Wiedemann-Franz law [79]. Therefore, a good thermoelectric material requires low thermal 
conductivity without compromising electrical conductivity. Ke is calculated using the BoltzTrap2 code 
as described in Eq. 38 in Chapter | in the framework of the two-current model [82]. Conversely, «1 is 
determined by the Slack equation. This equation specifically accounts for the heat conduction enabled 


by acoustic phonons through an anharmonic flip-scattering process and can be expressed as [84]: 


me} v3 


Kk, =A Tn 


(eq.40) 
Herein M is the average atomic mass, V is the average volume per atom, n is the number of atoms in 
the unit cell, and y is the Griineisen parameter calculated by the expression proposed by Julian [87] as 
9-12.24)? 
y= Tae? in which v; and v; are the longitudinal and transversal sound velocities, respectively. 
Vv 
The constant factor A is determined by the Griineisen parameter and can be calculated using the formula: 


2.43x107 
A= _ 0514, 0.228 - 


Figures 23(c-d) show the evolution of «ke and kl as a function of temperature for all series of 
materials discussed in this study. The curves representing ke in Figure 23(c) show a significant increase 
with increasing temperature for both p-type and n-type semiconductor materials in their pristine and 
HSC forms. This is due to the increase in electron energy, which increases the mobility of the charge 


carriers and consequently leads to higher electrical conductivity. Notably, the p-type alloys exhibit 
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higher ke values compared to their corresponding n-type counterparts. In contrast, for the HME alloy, 


ke decreases rapidly and remains constant at high temperatures. 


For comparison, p-type SrS: (Fe, K) has the highest xe, followed by SrS: (Fe, Li), then SrS: (Fe, 
Na), SrS and finally the SrS: Fe compound at 1200 K. The order is different for n-type carriers, where 
they are SrS: (Fe, K), SrS, SrS: (Fe, Li), SrS: (Fe, Na) and finally SrS: Fe. The behavior of «ke mirrors 
that of o and follows the proportionality presented in the Wiedemann-Franz law [79] for p-type carriers. 
For the n-type carriers, all compounds follow the law except SrS, which is consistent with Rosenberg's 
observation [88] that the Wiedemann-Franz law holds at both low and high temperatures, but may not 
hold at intermediate temperatures. Our results are not consistent with those of Hou et al. [80], which 


could also be due to the differences in the parameters used. 


To gain deeper insights into the temperature sensitivity of thermal conductivity, we plotted the 
curves of «1 in Figure 23(d), which clearly shows that «l decreases with increasing temperature. From 
100 K to 1200 K, the «l decreases from 126.51 W/m.K to 10.54 W/m.K for pristine SrS, from 69.04 
W/m.K to 5.76 W/m.K for SrS: (Fe, Li), from 60.92 W/m.K to 5.04 W/m K for SrS: (Fe, Na), from 
53.12 W/m K to 4.42 W/m K for SrS: (Fe, K ) and from 9.64 W/m K to 0.82 W/m K for SrS: Fe. This 
trend shows that higher temperatures lead to greater phonon scattering, resulting in reduced kl. This 
behavior of suppressed lattice thermal conductivity is desirable for improving the thermoelectric 


efficiency. 


Our original values are overestimated compared to previous literature work (2.60 W/m.K at 1200 
K) [80]. This can be attributed to the use of ShengBTE [87], which uses a fully iterative solution of the 
BTE and considers second and third order force constants instead of the Slack equation. This 
overestimation of «l can be clarified by considering the influence of inaccurately estimated y in the Slack 


equation, as mentioned in [88]. 


d. Figure of Merit 


Using the collected o, S, xe, and «1 data, we estimated the dimensionless quality factor (ZT) to 
evaluate the practical thermal efficiency of our compounds. The temperature-dependent ZT 
diagrams are shown in Figure 23(e). Upon careful analysis, we found that the ZT values for p-type 
SrS, SrS:Fe, and SrS:(Fe, p°) compounds outperform their n-type counterparts within the same 
material, reflecting the superior thermoelectric performance of p-type compounds underlines in 
comparison to their n-type carriers. Furthermore, it is evident that temperature has a significant 
influence on the thermoelectric behavior of these compounds. At 300 K, ZT remains relatively low 
and has values of 0.04, 0.24, 0.07, 0.02 and 0.3 for pristine SrS, SrS: Fe, SrS: (Fe, Li), SrS: (Fe, Na) 


and SrS: (Fe, K) in p-type doping carriers. However, when temperatures increase to 1200 K, the 
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maximum ZT values reach 0.24, 1.51, 3.11, 3.34, and 3.03 for p-doped carriers, maintaining the 
same order. In contrast, the maximum ZT value for n-doped carriers reaches 0.06, 0.02, 0.22, 0.08 
and 0.24 at 1200 K. This discrepancy can be due to the significantly higher electrical conductivity 


and electronic thermal conductivity of p-type carriers compared to their n-type counterparts. 


Furthermore, these impressive ZT values, which exceed one for p-type carriers, are largely due 
to the semiconductor nature of the narrow bandgap co-doped alloys. These results suggest that the 
HSC compounds are promising for high-temperature thermoelectric applications, especially 
compared to the HMF alloy. Furthermore, it shows that co-doping is more effective than mono- 
doping to achieve the best thermoelectric performance for binary SrS. It should be emphasized that 
the ZT values we calculated for SrS show slight deviations from those reported by Hou et al. [80] 


(0.08 for p-type and 0.15 for n-type). These deviations were explained in the respective section. 


3.2.4 Conclusion 


In the second section of this chapter, we examined an in-depth study of the three newly designed p° 


alkali metals (p°= Li, Na, and K) co-doped into SrS:Fe compounds at a selected doping level of 12.5% 
using rigorous ab-initio calculations within the PBE and PBE+mBJ methods. This study included a 
comprehensive investigation of the effects of co-doping on various properties, including structural, 
electronic, magnetic, optical, and thermoelectric properties. The key findings and conclusions from this 


endeavor are presented as follows: 


1. The results of the structural properties obtained with the PBE-GGA approximation have shown 
that the lattice constant a (A) increases gradually with the atomic number (Z), while the bulk 
modulus B (GPA) progressively decreases. Substitution of alkali metals in SrS:Fe has been 
shown to increase the overall stability of the system, facilitate the doping process, and strengthen 
ionic bonds. Calculations of the differences in total energies (AE) confirmed the stability of the 
co-doped systems in the ferromagnetic (FM) state, with a higher Curie temperature (Tc) being 
observed compared to room temperature. Among the co-doped systems, SrS: (Fe, K) in 


particular has the highest Tc, making it a promising candidate for spintronic applications. 


2. The electronic properties showed that all compounds were characterized as p-type 
semiconductors (HSCs), and it was observed that the energy gap became smaller in the co- 


doped systems compared to the semi-metallic (HM) SrS:Fe compound. 


3. In all compounds, we calculated a total magnetic moment of 5 Ls, which is increased compared 
to mono-doping, and a low local magnetic moment was created at the non-magnetic sites due to 


sp-d hybridization. 
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4. The analysis of optical properties in terms of dielectric function, absorption coefficient and 
optical conductivity has shown that all three types of alkali metal ions have the ability to enhance 
the infrared absorption of the SrS:Fe system and generate new peaks in the visible region and 


leads to a redshift. This property makes the materials ideal for solar cells applications. 


5. The analysis of the thermoelectric properties clearly showed that the p-type ZT values were 
significantly superior to those of the n-type ZTs. In particular, the p-type SrS: systems (Fe-Li, 
Na, and K) were found to be promising candidates for efficient thermoelectric materials, 


showing ZT values above 1 at 1200 K. 


In summary, Chapter 3 with its two parts has provided comprehensive insights into one of the two 
methods, namely co-doping. This method proves to be a powerful tool for tailoring the properties of 


HMFs to meet precise application requirements. 


In the next chapter of this dissertation, we will examine the second method, which focuses on 


dimensionality reduction. 
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Chapter 4 


Effect of Dimensionality Reduction on the Fe- 


Doped SrS Properties 


This part of the work is in progress. 


n the final chapter of this dissertation, we have studied the spin-resolved structural, electronic, and 
magnetic properties of Sr:.,.Fe,S monolayers at different concentrations (x = 0, 0.125, 0.25, 0.50, 
0.75, and 1). We have also compared the findings for these monolayers with those of their bulk 
counterparts, as examined in the initial section of Chapter 3, to get the effect of dimensionality reduction 


on the different properties. 


1.1 Introduction 


Two-dimensional (2D) materials have distinct properties that can differ significantly from their bulk 
counterparts due to their increased surface-to-volume ratio. This characteristic has provided a new 
foundation for hybrid device engineering, allowing us to explore and customize the materials’ 
exceptional properties. Although extensive research is being conducted on the applications of 2D 
materials, many investigations are still in their early stages, and there remain numerous unresolved 
issues that require attention. In the realm of 2D materials, one of the most notable examples is monolayer 
graphene [1], which captured the significant interest of material scientists since 2004, as we discussed 
earlier in Chapter 1. The SrS structure, akin to graphene, exhibits a large indirect bandgap in the near- 
visible spectrum and a remarkably flat band in the upper valence band [2]. This rendered it particularly 


intriguing in the post-graphene era for various applications. 


Keeping in view of these unique characteristics, in this chapter, we systematically investigate the 
structural, electronic and magnetic properties of SrS-to-FeS monolayers by varying the strontium ratio 
concentration using the powerful first-principle calculations. Until now, the study of Sri-Fe,S alloys 


has been limited to the bulk phase, as conducted for the first time in the previous chapter, with no prior 
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reports on the monolayer forms of Sri-xFe;xS alloys. Therefore, this study will forms a base for future 


theoretical and experimental investigations in this field. 
1.2 Computational Details 


The calculations presented in this chapter were performed using the Vienna Ab-initio Simulation 
Package (VASP) program [3, 4], within the Spin-Polarized Density Functional Theory (SPDFT) 
formalism and the Projector Augmented Wave (PAW) potential [5, 6]. The electron exchange- 
correlation functional was described by the Perdew-Burke-Ermzerhof (PBE) form of the generalized 
gradient approximation (GGA) [7]. To ensure accuracy, potential underestimations in bandgap and 
magnetic moment within GGA were cross-verified using the non-local Heyd-Scuseria-Ernzerhof hybrid 
functional (HSE06) [8]. In this approach, the electron-electron interaction was determined by a 


proportion of Hartree-Fock exchange with a = 0.25 and a default screening parameter of 0.2 A“. 


To create Sri-.Fe,xS alloys, the Sr atoms were systematically replaced with Fe atoms in varying 
ratios (x), corresponding to x = 0, 0.125, 0.25, 0.50, 0.75, and 1. For each alloy ratio, all possible 
combinations of Fe atom in the SrS host matrix were considered to select the most stable configuration. 
To model the different crystal structures of the systems, a (2 x 4 x 1) supercell of 16 atoms containing 
one Fe atom for x = 0.125 and a (2 x 2 x 1) supercell of 8 atoms containing one, two, and three Fe atoms 
for x = 0.25, 0.50, 0.75, respectively, were constructed based on a conventional SrS unit cell. Interlayer 
interactions between successive unit-cells were avoided by assuming a vacuum spacing of 20 A.A 
gamma scheme was used to sample the Brillouin zone (BZ), with a 12 x 12 x 1 k-point mesh in the unit 
cells and scaled accordingly to 6 x 3 x 1 and 6 x 6 x 1| k-point meshes, respectively, for the supercell 


simulations. 


The ground state was optimized with an energy convergence criterion of 10° eV between 
consecutive iterations. Additionally, the maximum Hellmann-Feynman force on each atom was 
constrained to be less than 103 eV/A during ion relaxation. To ensure accuracy, a cutoff value of 500 


eV for the plane-waves energy was employed. 


For charge distribution analysis, Bader analysis [9] was applied to different concentrations. The 
dynamic stabilities of bare SrS and FeS monolayers were assessed via phonon band dispersion using the 
PHONOPY code [10]. A (5 x 5 x 1) supercell was used for each monolayer with a 9 x 9 x | k-point 


mesh. 


Ab-initio molecular dynamics simulation (AIMD) [11] was employed to investigate the thermal 
stabilities of the monolayers. This was conducted using a supercell of (6 x 6 x 1) for A-SrS and (5 x 5 x 


1) for s-FeS at 300 K temperature, employing time steps of 1 fs. 
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1.3 Results and Discussion 
4.3.1 Structural Properties 


The SrS monolayer is an atomically thin material, composed of a single layer with a unique planar 
hexagonal honeycomb lattice structure. This lattice is entirely constructed from Sr-S bonds, exhibiting 
space group N°187 that falls within the D3n point group symmetry, akin to other well-studied monolayers 
like graphene, silicene, and BN. The unit cell of the SrS monolayer encompasses one Sr atom and one 


S atom. 


In contrast, the FeS monolayer possesses a distinctive square crystal structure within the P4/mmm 
space group N°123. Its unit cell comprises a total of four atoms, including two Fe atoms and two S 
atoms. Figure 1 offers an illustrative representation of the top and side views of both the A-SrS and s- 


FeS configurations. 


Figure 1. Top and side views of the repeated unit-cells of (a) the bare h-SrS and (b) s-FeS 


monolayers. 


The structural parameters including the lattice parameters a (A) and bond lengths d (A) of the 4-SrS 
and s-FeS are tabulated in Table 1. The lattice parameters were determined to be a = b = 4.837 A for h- 
SrS and 3.557 A for s-FeS. On the other hand, the bond lengths were determined to be ds,-s= 2.793 A 
for h-SrS, and 2.150 A for s-FeS. These values are in good agreement with those reported in the literature 
[12-15]. When comparing the lattice parameter of rock-salt SrS (6.059 A) to its hexagonal structure, we 
observed that it is larger in bulk form than in monolayer form. This can be attributed to dimensionality 
reduction, where interatomic interactions result in a different equilibrium spacing compared to a 


monolayer. 
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Table 1 


The equilibrium optimized structural parameters of bares h-SrS, s-FeS, and doped monolayers Sr1.xFe,S: 


lattice constant (a), shortest bond length (d), formation energy(E;), and elastic constants (Cj). 


Materials a=b(A) d(A) — Eyor(eV) = Cu (N/m) Cn (N/m) Cos (N/m) 
h-SrS 4.837 2.793 -1.402 28.993 23.801 - 
4.76* 2.75 - - = 2 
4.85° 2.80° - - = 2 
Sro.s7sFeo.125S 4.756 2.254 -1.176 4 s - 
Sro.7sFeo.25S 4.581 2.283 -1.049 - - - 
Sro.s0F eo.s0S 4168 2.217 -0.454 = 2 z 
Sro.2sFeo.7sS 4.085 2.232 -0.021 : : 2 
s-FeS 3.557 2.150 -0.427 T7575 46.595 60.542 
3.56° 2.16° - 7 7 : 
3.594 - - : - E 


“Ref [12], ° Ref [13], ° Ref [14], ‘ Ref [15]. 


We have evaluated the energy required for constructing the atomic assembly of both A-SrS and s- 


FeS monolayers by calculating the formation energy (E) using the formula as indicated in Chapter 3. 


The calculated EF 'f values are also indicated in Table 1. As it is seen from Table 1, the E 'f values for 
h-SrS and s-FeS are both negative, indicating that the bare monolayers are exothermic and unlikely to 
decompose after formation. Moreover, the formation energy of h-SrS (-1.402 eV) is notably less than 
that of s-FeS (-0.427 eV), indicating that h-SrS provides a more favorable environment for doping 
compared to s-FeS. Hence, when preparing Fe-doped SrS samples, it is advisable to create an h-SrS 


environment to ensure efficient and effective doping in practical experiments. 


Utilizing the Density Functional Perturbation Theory (DFPT) approach [16], we assessed the 
dynamic stability of the A-SrS and s-FeS monolayers by examining their phonon dispersions. In Figure 
2, the phonon spectra are presented along high-symmetry lines in the first Brillouin zone. It is evident 
that the phonon branches do not exhibit any imaginary frequencies, indicating that both structures 
remain dynamically stable under ambient conditions. However, slight occurrences of negative 


frequencies can be addressed by further increasing the supercell size. 


The phonon plots reveal that the h-SrS monolayer possesses six phonon branches, consisting of 
three acoustic and three optical branches. Conversely, the primitive cell of s-FeS, which contains four 


atoms, yields twelve phonon branches, comprising three acoustic and nine optical branches. The 
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maximum frequency in these monolayers reaches up to 140 cm! for h-SrS and 200 cm! for s-FeS, 
aligning with prior research findings [13, 14]. Our observations also indicate that the phonon spectra of 
both systems exhibit non-degeneracy along the directions of the Brillouin zone. Furthermore, we note 
an important phonon gap in the optical band, measuring approximately 50cm"! and 15cm‘ for h-SrS 


and s-FeS, respectively. 


Significantly, this phonon gap is most prominent in /-SrS and least in s-FeS, primarily attributed to 
disparities in the arrangement of atoms within their primitive cells. This phonon bandgap has the 
potential to enhance thermal conductivity through specific phonon-phonon scattering processes, 


suggesting that s-FeS may exhibit higher thermal conductivity compared to h-SrS. 


Subsequently, Ab-initio Molecular Dynamics (AIMD) simulations were conducted to perform a 
comprehensive evaluation of the thermal stability of the materials. Figure 3 depicts the energy 
fluctuations of the h-SrS and s-FeS bare monolayers over a 3000 fs period at 300 K. The AIMD 
simulations revealed consistent energy values for the bare monolayers throughout the entire simulation 
duration. The minor energy fluctuations arise from the rippling effect induced by the increase in 
temperature (the change in energy interval is minimal). This observation strongly supports the assertion 
that both the h-SrS and s-FeS monolayers possess stable structures and exhibit good thermal stability at 


room temperature. 
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Figure 2. Phonon dispersion spectra of the bare h-SrS and s-FeS, respectively. 
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Figure 3. The variation of the total energy with time during AIMD simulation at 300k over a 3000 fs. 


After confirming the thermodynamic, dynamic, and thermal stabilities of the pristine monolayers, we 


employed the widely recognized strain-energy approach [17, 18] to evaluate their mechanical stability. 


ai-ao 


By systematically varying the strain ratio (€ = , ai and ap indicate strained and strain-free 


lattice constants, respectively) within the range of -3% < ¢ < 3% with an incremental step size of 0.01, 
we acquired corresponding energy (E) values. The resulting E vs. € data was fitted using a quadratic 
polynomial equation. This fitting process allowed us to extract the in-plane stiffness using the following 


equation: 


Cj =~, (eq.1) 


In this context, Ag denotes the equilibrium area of the system. A system is deemed mechanically stable 


when it adheres to the minimum Born criterion for elastic stability in 2D materials [19]. 
For hexagonal and square crystals, the mechanical stability criterion are the following [20]: 
Cu>O;  Cir> |Cr2| for hexagonal structure (eq.2) 
Cir >0; Coo >O3 Cra > |Cr| for square structure (eq.3) 


Upon reviewing the elastic constants detailed in Table 1, we can affirm their compliance with the 
aforementioned stability criteria, confirming hence the mechanical stability of these monolayers. It is 


worth noting that experimental validation and further theoretical investigations are still required to 
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validate these findings. Specifically, the elastic constant Ci; measures a material's resistance to linear 
compression along the x direction [21]. Our computations revealed that the Ci; values for the pristine 
monolayers are markedly elevated in comparison to the other elastic constants with the s-FeS boasts the 
highest value of 77.575 N/m, followed by h-SrS. This suggests that both systems are highly resistant to 
uniaxial stress in the x direction, and s-FeS is particularly highly incompressible in this direction due to 


its exceptionally large Cy; value. 


Following our previous work on bulk structures, we conducted spin-polarized computations 
involving the introduction of Fe atoms into the pure A-SrS structure at four key concentrations: 12.5%, 
25%, 50%, and 75% Fe. In order to ensure energetically stable configurations for the resulting alloys, 
we explored various arrangements of the transition-metal atoms into the h-SrS monolayer for each alloy 
type, evaluating the total energy per atom. The most energetically stable configurations for each alloy 


are illustrated in Figure 4. 


Sro50F@o.50S Sto25Feo.758 


Figure 4. The ground-state configurations of predicted energetically stable Sr;-,Fe,S monolayers (Sr 


atom is designated by brown color, S by yellow, and Fe by grey). 


After geometric optimization, it was observed that all the alloyed ground-state structures maintained 
their original hexagonal configuration, albeit with slight distortions in the shape of their unit cells. Table 
1 furnishes detailed information on the lattice constants, formation energies, and bond lengths for the 
doped monolayer alloys. It is noteworthy that the lattice parameters consistently decrease with 
increasing Fe concentration in h-SrS. This reduction can be attributed to the disparity in ionic radii 
between Fe and Sr atoms. Additionally, they exhibit the same trend as their bulk counterparts, with a 
noticeable decrease in values from bulk to monolayer structure. This is also explained by the interatomic 


interactions resulting from the dimensionality reduction. 


Concerning the bond lengths, the average distance between Fe and the nearest S atoms in the doped 
monolayers ranges from dye-s = 2.217 to 2.283 A, depending on the specific doping concentration. The 


Fe-S bonds in the doped monolayers experience compression compared to Sr-S bonds before doping. 
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This compression can be attributed to the relatively small difference in ionic radii between Fe and Sr 
atoms. However, it is noteworthy that the Fe-S bond length shows minimal variation with changes in 


the doping concentration, indicating the structural stability of the alloys. 


The formation energy for the doped monolayers has also been computed utilizing the previously 
outlined formula from Chapter 3. The outcomes, detailed in Table 1, demonstrate that Fe atoms can be 
effectively incorporated into the monolayer h-SrS through an exothermic process. Upon examining the 
formation energies, a gradual increase is observed as the doping concentration escalates from 0.125 to 
0.75. This pattern signifies that the h-SrS with 12.5% Fe doping displays the highest thermodynamic 


stability among the alloys, suggesting its most energetically favorable configuration. 


Upon comparing the formation energies of monolayers to their bulk counterparts (please refer to 
Table 1 in Chapter3), it becomes evident that the formation energies in 2D configurations are greater 
than those in 3D configurations. This leads to the conclusion that in a bulk structure, atoms establish 
more connections and form a higher number of bonds, ultimately resulting in increased binding energy 
and, consequently, lower formation energy. Moreover, the higher surface of 2D materials enables 
increased surface energy due to the incomplete coordination compared to the bulk counterparts; thereby 


contribute to a higher energy of formation. 


Also in terms of comparison, the stability demonstrated through co-doping SrS: Fe with alkali 
metals surpasses that achieved by reducing its dimensionality. This highlights the pronounced 


effectiveness of co-dopage method in enhancing stability in our specific case. 
4.3.2 Electronic Properties and Bader Charge Analysis 
4.3.2.1 Band Structures Analysis 


In Figure 5, we have depicted the band structures of both the pristine compounds and Sri-xFexS 
monolayer alloys using the Heyd—Scuseria—Ernzerhof (HSE06) method. Additionally, we have 
tabulated the computed bandgap values for these compounds in Table 2, comparing them with those 
obtained via the GGA method. It is crucial to note that the GGA method, employed in our computational 
analysis, tends to underestimate the bandgap values of the materials under study. This tendency can 
potentially affect the precision of our findings, especially in characterizing the electronic and optical 
properties of the materials. To tackle this challenge, the HSE has emerged as a promising alternative 


due to its capability to offer more dependable and precise predictions of bandgap values. 


From Figure 5, it is evident that the A-SrS exhibits the characteristics of a non-magnetic 
semiconductor, displaying an indirect bandgap of 2.825 eV within the GGA method and 3.861 eV within 
HSE, which aligns with earlier researches [12, 13]. Notably, the CBM is positioned at the I’ point of the 
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BZ, while the VBM is located at the M point. It is obvious when the dimensionality of SrS decreases 
from 3D to 2D; the bandgaps maintain indirect with the gap values increase from 3.435 eV using the 
mBJ to 3.861 eV using the HSE for 2D monolayer. This is expected as the lattice parameter demonstrates 


the opposite trend when transitioning from 3D to 2D. 


In contrast, the electronic behavior of the s-FeS monolayer differs from that of h-SrS. The analysis 
reveals a metallic character, with the Fermi level intersecting both the valence and conduction bands, 
resulting in a non-magnetic metallic state according to both the GGA and HSE methods. These 
observations validate a previous study [14], underscoring the reliability and accuracy of the obtained 


results. 


The introduction of Fe atoms through doping in the pristine h-SrS monolayer induces significant 
changes in its band structure. This doping introduces impurity levels within the initial bandgap, causing 
a notable shift in the electronic character towards a HSC character within the HSE. Notably, this 


alteration induced by doping leads to a distinct lowering of the position of the CBM. 


As detailed in Table 2, within the GGA method, the bandgap diminishes from 2.825 eV in its pristine 
form to a range of 1.822—0.000 eV in up-spin channel and to 0.664-0.000 eV in down-spin channel, 
depending on the doping concentration of the Fe atoms. Similarly, the HSE method predicts a decrease 
in the bandgap from 3.861 eV to a range of 3.225—3.049 eV in up-spin channel and to 2.935-2.297 eV 


in down-spin channel. 


The PBE results indicate that the HSC monolayers including, Sro7sFeo.125S, Sro.7sFeo25S, and 
Sro.soFeo.soS, display total indirect bandgaps in [H-I'], [K-I'], and [I-K], respectively, while the 


Sro.25Feo.75S demonstrates metallic character. 


However, under the application of the HSE method, it is observed that the Sro.2sFeo7sS transition 
from a metal to HSC with a total direct bandgap of 2.263 eV in [I- I] direction, while the remaining 


alloys preserve their semiconducting properties within the same bandgap direction. 


It is important to highlight that the effect of doping on the band structure, as observed with both the 
PBE and HSE methods, shows non-linear behavior with respect to the doping concentration. The 
variation of the total band gaps, defined as the smallest distance between the CBM and VBM of both 


spin channels, is tabulated in Table 2 and shown in Figure 6 for visualization. 


The significant reduction in the bandgap of h-SrS after the incorporation of Fe atoms can be 
attributed to the hybridization between the orbitals of the Fe dopant atoms and the neighboring Sr/S 


atoms, as further explained by the partial density of states (PDOS) analysis in the following section. 
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Comparing the results for doped monolayers with their bulk counterparts, we conclude that the 
reduction in dimensionality shifts the nature from HMF to HSC by creating a gap in the spin-down 
channel. In addition, the doped monolayers exhibit a particularly flat valence band, a typical feature of 
their 2D configurations. This flat valence band represents an important and advantageous feature to 


consider in thermoelectric applications. 


4.3.2.2 Bader Charge Analysis 


Bader charge analysis provides valuable insights into the redistribution of electron charges within 
a material. Here, we performed Bader charge analysis to understand the partitioning of the total charge 
of the h-SrS monolayer into localized atomic charges after Fe dopant substitution. The results are 


summarized in Table 2. 


The observed variation in charge transfer values is consistent with the corresponding trend of 
bandgap variation within the monolayers calculated using the GGA method. In the original h-SrS 
monolayer, both Sr and S atoms carry a net Bader charge of 1.236e, reflecting the electron transfer from 
Sr to S atoms caused by their different electronegativities. When doping with Fe atoms, the amount of 
charge transferred to the Fe atom varies from 0.399e to 0.575e to 0.609e to 0.268e for Sro.g7sFeo.125S, 


Sro.75Feo.25S, Sto.soFeo.soS, and Sro.2sFeo.75S, respectively. 


These values are lower than those in bare h-SrS, indicating a less robust interaction between Fe and 
the neighboring S atoms compared to that between Sr and S atoms. A reduced Bader interaction implies 
a weaker bond between the doped Fe ions and the neighboring atoms. This could possibly be due to 


changes in the electronic structure or atomic arrangement caused by the doping process. 
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Figure 5. Spin-polarized electronic band structures of the bare h-SrS, s-FeS and Sr-x FexS monolayer 


alloys calculated with the HSE06 method. The Ef is set to zero. 
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Figure 6. Visualization of the bandgap fluctuation of the bare h-SrS, s-FeS and Sr-x FexS monolayer 


alloys. 
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Table 2 


The calculated electronic properties of bare h-SrS, s-FeS, and doped monolayer alloys Srj-xFe;S, 
including bandgap values for up-spin channel, down-spin channel, and total gap within both PBE and 


HSE functionals, along with Bader charge (Aq) transfer from Sr (Fe) to S atoms. 


Materials Up-spin channel Dn-spin channel Total gap Aq 
direction Value(eV) direction Value(eV) direction Value(eV) (e) 
PBE [MT] 2.825 [M-I’] 2,825 1.236 
[M-I'}* 217" 
[M-I']> 2.54? 
HSE [MI] 3.860 [M-I] 3.860 
[M-T] 3.74? 
Sros7sFeo.12s8 PBE [X-I] 1.822 [H-T] 0.229 [H-I] 0,229 0.399 
HSE [X-T] 3.154 [H-T] 2.451 [H-T] 2.435 
Sro7sFeo2s8 PBE [M-I'] 1.803 [K-T] 0.664 [K-I] 0.664 0.575 
HSE [M-T] 3.105 [K-I] 2.935 [K-I] 2.830 
SrosoFeosoS PBE [I-I] 1.780 [I'-K] 0.198 [T-K] 0.198 0.609 
HSE [I -T] 3.050 [I-K] 2.564 [I-K] 2.564 
Sro2sFeo.758 PBE Metal Metal Metal 0.268 
HSE [IT] 3.225 (T-T] 2.297 [T -T] 2.263 
PBE Metal Metal 0.771 
Metal‘ 
HSE Metal Metal 


*Ref [12], °Ref [13], ‘Ref [14]. 


4.3.2.3 Total and Partial Densities of States Analysis 


The TDOS and PDOS of the bare 4-SrS and s-FeS are visually depicted in Figure 7. 

The analysis of the TDOS for the h-SrS bare monolayer clearly affirms its semiconducting nature. 
This is evident from the absence of energy states at the Er level, signifying the presence of a bandgap 
between the CB and VB. Moreover, the examination reveals a distinct two peaks precisely situated at 
the VB edge at -0.117 eV and -0.932 eV, respectively, representing the highest occupied energy levels 
in the material. Further investigation via the PDOS clarifies that the VB arises mainly from the 3p, 
orbitals of S atom and a weak hybridization between the Sr: 5s-orbitals and the S: 3p;-orbitals, with the 


S orbitals providing the most substantial contributions. Conversely, the emergence of the CB in the 
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energy range of 3.860 to 6 eV is predominantly attributed to the 5s-orbitals of the Sr atoms, accompanied 


by a minor contribution of the S: 3p,-orbitals. These observations align with prior studies [13]. 


The TDOS of the s-FeS bare monolayer displays a different characteristic from that of h-SrS, within 
a sharp peak extending towards the CB edge giving a metallic character. Analysis of the PDOS of s-FeS 
reveals that the VB in the energy range from -6 eV to 4 eV is dominated by both the S: p-orbitals and 


the Fe: d-orbitals, within the highest contribution coming from the S: 3p, and Fe: 3dx2-y2-states 


In the energy range from -3.5 eV to 2.5 eV, the contribution of the S: p-states decreases and Fe: d- 
states dominate. Specifically, Fe: 3dxy dominates in the energy range of -3.5 eV to -2 eV, while Fe: 3dz” 
and Fe: 3dx?-y” dominate in the energy range of -2 eV to 1 eV, intersecting the Er level. This confirms 
the metallic characteristic exhibited by this bare monolayer. From 1 eV to 2.5 eV, both Fe: 3dx, and Fe: 
3dz° are dominant. A minor contribution in this range comes from the S: 3p,. These findings are in 


agreement with the results presented in reference [14]. 
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Figure 7. Electronic total (TDOS) and partial (PDOS) densities of states of the bare h-SrS and s-FeS 


monolayers. 


The TDOS and PDOS of the ternary-doped monolayers at different concentrations are displayed in 


Figures 8 and 9, respectively. 


Looking carefully at these numbers, the disparity between spin-up and spin-down TDOS in the Sri- 
xFex,S monolayers is striking. This asymmetry is in perfect agreement with the ferromagnetic ground 
states of all systems and confirms their magnetic nature. Furthermore, the presence of the impurity states 
in the minority VB, located in close proximity to the Er, confirms their p-type character even after 
doping with impurities, which is fundamental for enabling a wide range of electronic devices and 


technologies. 
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Specifically, in the majority-spin, the Fe: 3d states and S: 3p states of the four compounds hybridize 
with each other. This leads to the preservation and reduction of the bandgap compared to that of the bare 
h-SrS. In the minority-spin, the p-d hybridization also shifts the 3d-Fe and 3p-S states towards the Ep, 
while maintaining the bandgap. This results in a distinctly HSC character for these compounds. This 
behavior is definitely different from that of their bulk counterparts, where the p-d hybridization and the 


emergence of 3d states in the minority states at Er lead to a HMF character with 100% spin-polarization. 


We come to a very important conclusion in our case that is by reducing the symmetry in 2D 


monolayer form we can drastically change the character of the compounds. 


It is also evident from the Figure 8 that with increasing the concentration of the Fe impurities from 
12.5% to 75%, respectively, the p-d hybridization become notably pronounced and the 5s states effect 


decreases, which aligns also with the bulk results. 
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Figure 8. Electronic TDOS of Srj-.FexS (x= 0.125, 0.25, 0.50, and 0.75) monolayer alloys. The 


positive and negative values represent the spin-up and spin-down channels, respectively. 
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We further analyzed the PDOS of the monolayer alloys and plotted them in Figure 9. The analysis 
of the curves revealed that, for the Sro.s7sFeo.125S compound, the majority-spin VB, in the energy range 
of ~-5 eV to ~-1 eV is primarily derived from the 3px and 3p, states of S atoms, and the 3dxz, 3dyz, and 
3dz’ of the Fe atoms, with a minor contribution coming from the 5s-states of the Sr atoms. As the Fe 
content reaches 25%, 50%, and 75 % Fe, the 3d,, and 3dx, of Fe atom increase and become dominant 
compared to the 3px and 3p, states of S atoms and 5s-states of the Sr atoms, which decrease with 


increasing the concentration. 


In the majority-spin CB, in the energy range of ~-3 eV to around ~6 eV, the Sr: 5s become dominant 


through all the four concentrations, with a small contributions coming from 3p, and 3p; of the S atoms. 


For the minority-spin channel, the VB is mainly dominated by the S: 3p, and S: 3p, states and the 
VBM at around -1 eV is dominated by Fe: 3d,, and Fe: 3d,, states for 12.5% Fe, by the Fe: 3d, states 
for 25% Fe, the Fe: 3d,2 states for 50% Fe, and by the Fe: 3d,2 and Fe: 3dy, for 75% Fe. 


The CB on the other side, in the energy range of ~-3 eV to ~5 eV, is dominated by the Fe: 3dx, 


states and the Fe: 3dx?-y’ and a little contribution coming from the 5s states of the Sr. 


Overall, the impact of the p-d hybridization on the electronic behavior is even seen at low- 
dimensionality, underscoring the significance of this phenomenon in influencing the bandgap and 


changing the character of materials. 
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Figure 9. Electronic PDOS of Srj-,FexS (x= 0.125, 0.25, 0.50, and 0.75) monolayer alloys. 


4.3.3. Magnetic Properties 


The magnetic properties of the Fe-doped h-SrS monolayer alloys are evaluated based on their 


magnetic moments, and the outcomes are presented in Table 3. 


The results revealed an integer value of the magnetic moment (Mror) of 4g per supercell across 
all concentrations. These findings align with the bulk structure results reported earlier in Table 5 of 
Chapter 3. The distribution of magnetic moments within the Sr;-,FexS monolayers conforms to Hund’s 
rule and crystal field theory. Interaction with neighboring S ions induces a crystal field splitting, causing 
the 3d- orbitals of the transition metal atom Fe to separate into distinct energy levels. In this context, the 
trigonal-bipyramidal arrangement of Fe introduces a symmetry that can be depicted using the D3n point 
group. These levels encompass a non-degenerate orbital 'a' (dz’), a two-fold degenerate orbital 'eg' (dxy 


and dx2-y2), and another two-fold degenerate orbital 'eg»' (dx, and dy,), as illustrated in Figure 10. 


In the spin-up configuration of Fe** (d°), the ‘a’, 'egi', and 'egs' orbitals are situated below the Er 
level and are fully occupied, while the 'a' orbital in the spin-down channel lies above the Er. This 


electronic configuration can be represented as a t | ego Tt egi Tt, resulting in a Mror of 4p. 
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Figure 10. Schematic diagram representation of the crystal field splitting of the Fe** sites. Modified 
from [22]. 


Upon comparing the Mror per supercell with the local magnetic moments per atom (Table 3), it 
becomes evident that the primary contribution to the magnetic moment arises from the Fe atom, 
affirming the findings presented in the PDOS analysis. Smaller magnetic moments are attributed to the 
nearest neighboring Sr and S atoms, which exhibit ferromagnetic coupling with their adjacent Fe atoms 


due to their positive signs. 


Furthermore, it was observed that the local magnetic moment of Fe displays a linear variation, where 
it increases with increasing doping concentration from 12.5% to 75%. This behavior suggests a strong 
interaction between the Fe dopants and the adjacent surrounding atoms, reinforcing the ferromagnetic 
order in these monolayers. It should be mentioned that this trend differs from that of the bulk structures, 
where the M™ shown a decrease with concentration increasing. This may be also attributed to the 
difference in the arrangement of the atoms and how they interact in 2D monolayer structure, where there 


are fewer neighboring atoms, compared to the bulk structure. 


Table 3 


Total magnetic moment (Mror) and local magnetic moment on each Sr, S, and Fe atoms of Srj-xFexS 


(x= 0.125, 0.25, 0.50, and 0.75) monolayer alloys. 


Concentration (x) M“ (us)  MF(us) MS*(up) = MS(up) 


0.125 4.000 3.545 0.005 0.028 
0.250 4.000 3.595 0.011 0.035 
0.500 4.000 3.556 0.028 0.051 
0.750 4.000 3.587 0.089 0.037 
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1.4 Conclusion 


In this chapter, we have done a spin-resolved ab-initio study of the 2D monolayer form of Sri- Fe,S 
(x = 0, 0.125, 0.25, 0.50, 0.75, and 1) ordered alloys using a plane-wave pseudopotential DFT method 
within the VASP code. This study involved a comprehensive investigation of the effect of 
dimensionality reduction on different properties including the structural, electronic, and magnetic 


properties. The principle outcomes and conclusions are outlined as follows: 


1. The outcomes of formation energies, molecular dynamics simulations, and elastic stiffness 
constants indicated that each of the examined bare materials, SrS (FeS), is energetically 
favorable and remained thermally and mechanically stable under ambient conditions. Phonon 
spectra calculations further demonstrated that both materials exhibited dynamic stability in their 
lattice structures. In terms of ground state results, it is found that SrS maintained stability in a 
hexagonal structure akin to graphene, while FeS shown a preference for a square structure, 
aligning with prior theoretical studies. When Fe impurities were introduced in the Sr-sites 
through doping, all four monolayers retained their structural and thermodynamic stability in the 
hexagonal form, as confirmed by the formation energy. Upon comparing the formation energies, 
it is determined that the h-SrS system with 12.5% Fe doping displayed the lowest energy among 


the doped cases, indicating its most stable composition. 


2. The electronic properties analysis confirmed that both the PBE and HSE06 methods identified 
the bare monolayer h-SrS as a non-magnetic semiconductor with an indirect bandgap. On the 
other hand, the s-FeS system was identified as a non-magnetic metal. Concerning compositions, 
according to PBE results, three compounds namely, Sro.7sFeo.125S, Sto.75 Feo.25S and Sro.soFeo.s0S 
exhibited the behavior of HSCs, while Sro.2sFeo.7sS behaved as a metal. However, when the 
HSE06 method was applied, Sro.25Feo.75S transformed into a HSC with a direct bandgap in both 
spin-directions, while the other alloys maintained their semiconducting properties with an 
enhancement in bandgap values. We have concluded that the dimensionality affected the 


electronic character of these compounds that shown the HMF in bulk form. 


3. In all doped monolayers, we calculated a total magnetic moment of 4 up, which is consistent 
with the bulk structures. Additionally, low local magnetic moments on the non-magnetic sites 
were generated due to the p-d hybridization. We also observed that the local magnetic moments 
on Fe sites increased with increasing concentration, contrary to their bulk counterparts, 


confirming the strong ferromagnetic character at low dimension. 
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General Conclusion and Future Perspectives 


n the pursuit of new materials aimed at optimizing energy usage, advancing information storage, 
and promoting sustainability, spintronics, optoelectronics, and thermoelectrics have captured 
significant attention from researchers, being regarded as pivotal technologies for the future. 
Identifying materials with the requisite properties for effective application across these three domains 
proves to be a challenging task, owing to the specific criteria set by material properties. Presently, there 
exists a pressing need to develop new materials that are both more abundant, less harmful, and more 
cost-effective than their current counterparts. In light of this, the research approach has focused on 
enhancing the different properties of materials that already exhibit good ferromagnetic properties for 


spintronics. 


Dilute magnetic semiconductors (DMSs), the building blocks of spintronics, can be obtained by 
doping transition metals into semiconductor host matrices. Due to their high Curie temperature, low 
cost, ease of synthesis, and high electrical conductivity, iron (Fe) transition metal atoms are a great 


option for ferromagnetic dopants used in this dissertation to improve host matrix properties. 


Advances in theories and improvement of computational methods such as density functional theory 
(DFT) in recent years have led to significant advances in the understanding and analysis of the physical 
properties of materials. The work of this dissertation relied on theoretical research using ab-initio 


calculations to contribute to the study of new materials desirable for the above applications. 


First, we investigated the structural, mechanical, electronic and magnetic properties of Sri-x FexS 
compounds (x = 0, 0.125, 0.25, 0.50 and 0.75) using the FP-LAPW method and the implemented PBE/ 
PBE+mBJ approximations in the WIEN2K code. Our aim was to understand the influence of Fe doping 
on various properties of binary SrS and to explore new potential semi-metallic ferromagnets (HMF) for 
spintronic applications. Examination of the structural properties revealed that these materials were more 
stable in the ferromagnetic (FM) phase than in the anti-ferromagnetic (AFM) and non-magnetic (NM) 
phases. The formation energies indicated their thermodynamic stability. An increase in Fe concentration 
led to a decrease in the lattice parameter a (A), accompanied by an increase in the hardness B (GPA). 
Analysis of their mechanical properties showed that they were mechanically stable and showed greater 
resistance to unidirectional compression than to shear deformation. They showed high anisotropy and 


were stiffer in the [100] direction. They can be classified as ductile compounds. They also exhibited 
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ionic bonding nature and resist high temperature treatment. For the selected concentrations, the DMSs 
Sri-xFexS (xk = 0.125, 0.25 and 0.50) showed a total magnetic moment of 4 pg and a ferromagnetic half- 
metallic character (HMF) due to the strong p-d hybridization, making them good candidates for 
spintronics. In contrast, Sri-sFexS (x=0.75) showed a metallic character. Based on the obtained negative 
and positive exchange constants Noa and Nop, it was observed that the exchange coupling between the 
valence band of SrS and Fe-3d states is anti-ferromagnetic, whereas between the conduction band of 


SrS and Fe-3d states, the exchange coupling is a mixture of ferromagnetic and anti-ferromagnetic. 


The next step was to enhance the properties of these materials to make them more efficient. Co- 


doping and dimensionality reduction were the two methods used to accomplish this. 


The study of co-doping initially involved selecting the dopant concentration at x = 0.125 and 
incorporating the alkali metals Li, Na, and K alongside Fe. This was consistently done using the PBE 
/PBE+mBJ as incorporated in the WIEN2K code, complemented by the Boltzmann transport equations 
integrated into the BoltzTrap code. The investigation encompassed structural, electronic, and magnetic 
properties of the quaternaries with the formula Sro.37sFeo.o6p°o.06S (p° = Li, Na, and K), followed by an 
examination of their optical and thermoelectric properties. The outcomes regarding structural properties 
revealed a gradual increase in lattice constant with atomic number (Z), coupled with a progressive 
decrease in bulk modulus. Introducing alkali metal in SrS: Fe proved instrumental in maximizing the 
system's overall stability and reinforcing ionic bonds. Calculations of total energy differences (AE) 
affirmed the stability of the systems in the ferromagnetic state, with a notably higher Curie temperature 
(T,) observed compared to room temperature (RM). Notably, SrS: (Fe, K) exhibited the highest T, 
among these co-doped systems. Examination of the electronic and magnetic properties revealed that all 
materials exhibited p-type semiconductor behavior (HSCs). In particular, it was observed that the energy 
gap in the co-doped systems was smaller compared to the SrS:Fe compound. A total magnetic moment 
of 5 us was calculated, which represents an increase compared to mono-doping, and a small local 
magnetic moment was generated at the non-magnetic sites due to the sp-d hybridization. Examination 
of the optical properties, including real and imaginary parts of the dielectric function, absorption 
coefficient, and optical conductivity, revealed that each of the three alkali metal ions contributed to 
enhanced infrared absorption in the SrS: Fe system and the formation of new peaks in the visible 
spectrum. Such a property significantly increases the suitability of these materials for applications in 
solar cells. Finally, the analysis of thermoelectric properties, in terms of Seebeck coefficient, electrical 
conductivity, thermal conductivity and merit factor, highlighted the exceptional performance of the p- 
type ZT values compared to their n-type and mono-doped counterparts and exceeded them 1 at 1200 K. 


This makes them promising candidates for high-temperature thermoelectric applications. 


The dimensionality reduction method involved reducing the dimension of the ternary compounds 


Sti-x FexS (x = 0, 0.125, 0.25, 0.50, 0.75, and 1) from their 3D bulk structures to their 2D monolayer 
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counterparts using a plane-wave pseudopotential (PAW) method and the PBE/PBE+HSE 
approximations within the VASP code. This investigation included a thorough examination of their 
structural, electronic, and magnetic properties. The findings from the analysis of formation energies, 
molecular dynamics simulations, and elastic constants indicated that each of the examined pristine 
materials, SrS (FeS), had favorable formation energies and maintained thermal and mechanical 
stabilities under normal conditions. Calculations of phonon spectra also confirmed the dynamic stability 
of their lattice structures. Regarding their ground state configurations, SrS was observed to remain stable 
in a hexagonal structure like graphene, while FeS preferred a square structure, which is consistent with 
previous theoretical studies. The influence of dimensionality reduction is clearly observed in the 
electronic properties of the doped monolayers, which exhibited a behavior that differs significantly from 
the HMF character observed in the bulk counterparts. According to the PBE results, Sro.g7sFeo.125S, Sto.75 
Feo.25S and Sro.soFeo.s0S showed HSC behavior, while Sro.2sFeo.75S behaved like a metal. However, with 
the application of the HSE method, Sro2sFeo7sS transformed into an HSC with a direct bandgap in both 
spin directions. Finally, we calculated a total magnetic moment of 4 [1p in all doped monolayers, a result 
from p-d hybridization that is consistent with results in bulk structures. We also found an increase in 
local magnetic moments at Fe sites with concentration in contrast to their bulk counterparts, confirming 


the robust ferromagnetic nature at small dimensions. 


In summary, both co-doping and dimensionality reduction show promising key performances for the 


desired applications. 


In perspectives, we plan to calculate the optical and thermoelectric properties of our monolayers to 
enable a comprehensive comparison between 2D and 3D scales. Furthermore, it would be interesting to 
evaluate the reliability of other methods discussed in Chapter 1, such as inducing defects or external 


stimuli in SrS-based DMS. 
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Appendix. A 


This appendix contains basic equations that serve as initial steps in understanding the process of 


finding an appropriate basis set. 


A.1 Bloch’s Theorem 


In the field of solid-state systems, it is common for the units of interest to have specific translational 
symmetries in space. This property is suitable for the application of periodic boundary conditions (PBC) 
[1-3], a technique that significantly streamlines calculations for real systems. Under PBC, the system is 
encapsulated in a unit cell, systematically replicated throughout space and defined by the vectors 
Q@y, Az, and a3 [4]. This well-structured setup ensures that each grid point R can be precisely 


constructed, facilitating detailed analysis and calculations. 
R= nay NzAz + N303 (eq.1) 


The Bravais lattice is characterized by three integer coefficients, namely nj, nz and n3. This grid is 


paired with a reciprocal grid : 
G = m,b, mb, + mb; (eq.2) 


Each reciprocal lattice vector G can be expressed as a linear combination of the reciprocal basis 


vectors by, bz, and b3 where mj, m2, and m3 are integer values. These reciprocal basis vectors are 


derived from their counterparts in real space. 


ead 2 —> —> 
b, = = (az x a3) (eq.3) 
so 2 — —> 
b2 = = (is x ay) (eq.4) 
ad 2 — —» 
b3 = — (ay x az) (eq.5) 
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Here, 2 = ay(@z X @3) represents the volume of the unit cell in real space. Both the real and reciprocal 
basis vectors adhere to the relationship a,.b, = 26;;. The Schrédinger equation governing a single 


electron under periodic boundary conditions (PBC) is provided by [4]: 
1 = => = 
(- Ayz 4 v@)) WF) = ew?) (eq.6) 


Within a system governed by a periodic potential V(¥) that adheres to V(7' + R) = V(f), the Bloch 


theorem, credited to Swiss physicist Felix Bloch in 1929 [5], postulates that all eigenstates (7) of the 


one-electron Hamiltonian can be expressed as a composite of a plane wave phase factor e'" and a 


periodic function u, (7). This function mirrors the precise periodicity exhibited by the Bravais lattice 


[4]. 


Wn) = eu, (P) (eq.7) 
Ung? + R) = Unp(F) (eq.8) 


In this context, the quantum numbers 'n' and ‘k' refer to the band index and the Bloch wave vector in 
reciprocal space, respectively. It is a fundamental principle that any periodic function can be expressed 


as a sum of plane waves in a basis set [2]: 


Un (TY) = Da CnhkG elGr (eq.9) 


Consequently, the one-electron wave function can be expanded as follows: 


Wak’) = Xe Cake gibt (eq.10) 


In practical applications, the plane waves are typically truncated within a cutoff |Gnqx|, which 
corresponds to an energy cutoff Ecyt = |K + Gmax|*/2. Furthermore, the Bloch wave vectors 'k' can be 
confined to the Brillouin zone (BZ), which represents the first Wigner-Seitz cell in reciprocal space. 
This restriction is possible because all other Bloch states ‘k' can be calculated by combining a Bloch 


state within the BZ (Kgz) with an additional reciprocal lattice vector 'G' [4]. 


k= Kpz +G (eq.11) 


This implies that the Bloch states (eigenvalues and eigenfunctions) exhibit periodicity in reciprocal 
space [4]. To streamline computational efforts, symmetry considerations come into play. Specifically, 


only k-points within the irreducible wedge of the Brillouin Zone (IBZ) are taken into account, each with 
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a corresponding weight denoted as w,. With this approach, the summation over k-points for a periodic 


function F(k) spanning the entire Brillouin Zone can be simplified as follows: 
1 : 0 . 
ae Spy Ak F(k) = a Spy dk F(R) = DR? F(R) = Di? KF (k) (eq.12) 


Here, the k point grids and their associated weights w,, are typically generated using the Monkhorst and 
Pack scheme [6]. Qgz represents the volume of the Brillouin Zone, while Q denotes the volume of the 


real-space unit cell. 


The band structure of solids is represented by €,(k) a function dependent on the k vector for a 
specific band 'n'. The Density of States (DOS) for each band 'n' is determined by integrating over the 
IBZ [4]: 


Qf IBZ 
D,(©) = ora lay dk d(€ 7 En (k)) = dk w(€ ae €,(k)) (eq.13) 
As a result, the total DOS, encompassing contributions from all bands, is calculated as: 


D(€) = Yn Dr) (eq.14) 


A.2 Fourier Transformation 


In physics, the Fourier transform (FT) is a powerful mathematical operation that transforms a 
function and provides a detailed representation of the frequencies contained in the original data. Given 
a three-dimensional rectangular box (a finite system) with side lengths L,, Ly, and L,, any sufficiently 


smooth function f (7) that satisfies the periodic boundary conditions can be considered [7]. 
fF + Lye.) = f(F + Lyey) = f+ L,e,) =f) (eq.15) 
Eq.15 can be represented as a series of Fourier as follows: 
ee ik? 
f® = ork e (eq.16) 
Where k represents a wave vector within the reciprocal space of the system, which has a volume Q = 
LxLyLz. 


210 
+n 


ae ey, {ny, Ny, n,}e integer (eq.17) 


21 21 
k= My Tex Hay ey 


Subsequently, the Fourier coefficients F(k) can be computed through the following process: 
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F(k) = fidrf@e™* (eq.18) 
Some frequently utilized formulas for simplifying equations: 
ale ett = 5(7) 
far et" = 06(k) (eq.19) 
The principles and formulas previously discussed for the three-dimensional system can be extended 


and applied to a two-dimensional environment, taking the respective dimensions and coordinates into 


account. 
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Appendix. B 


This appendix provides a guide to better understand the WIEN2K and VASP codes and includes 
the flowchart of basic programs that serve as first steps in understanding the process of how the programs 


are used and executed. 


B.1 WIEN2K Code 


1.1 “Master Input™ File case.struct: Preliminary Requirement in WIEN2K 


To start DFT calculations in WIEN2K for a specific system, it is important to create a .struct file 
that represents the atom types involved in the simulation. This file, called “case.struct” (where “case” 
is the file name), serves as the main requirement in WIEN2K. The generation process is facilitated by 
StructGen, which is accessible via a web browser and the w2web interface or alternatively via the 
command line of an xterm. This crucial first step is to create a blank structure template into which we 
can enter relevant structure details such as atom symbols, lattice parameters, space group, atom 
positions, angles, atomic number (Z), and Rr values, which can be set either manually or automatically. 
After completing the StructGen process, the file "case.struct" is generated, which serves as the master 
input file for all subsequent programs. This step also automatically creates an input file with the atomic 


configurations, labeled “case. inst.” 


1.2 Main Programs in WIEN2K 


Once the two basic input files have been generated, the calculations are initialized using the “init- 
lapw” command. This command provides a user-friendly step-by-step guide to setting up the 
calculations and generating inputs for the main programs. By initializing the calculations, several 
automated steps are carried out seamlessly. Where 'x' represents the script for starting WIEN2k 


programs. 


e  xnn_This command calculates the nearest neighbors within a specified distance. 

e xsgroup_ This command calculates both the point and space groups for a given structure. 

e x-symmetry_ This command generates space group symmetry operations, determines the point 
group for individual atom locations, and generates the local moment extension (LM) for lattice 


harmonics and local rotation matrices. 


192 


Appendices 


e xIstart_ This command generates atomic densities and defines how the orbitals are treated in 
band structure calculations, including considerations such as core or band states, local orbitals, 
etc. 

e xkgen_ This command creates a k-network within the Brillouin zone (BZ). 


e xdstart_ This command generates an initial density for the Self-Consistent Field (SCF) cycle. 


After initializing the calculation, the SCF process is initiated and repeated until solution 
convergence is achieved. For non-magnetic materials, such as SrS in our case, the SCF cycle can be 
triggered with the “run_lapw” command. For materials with ferromagnetic properties, such as Fe-doped 


SrS in our scenario, the “runsp_lapw ” command is used. The SCF cycle includes the following phases: 


e LAPW0 (POTENTIAL): This phase generates the potential based on the calculated density. 

e LAPWI1 (BANDS): Here the valence bands, including eigenvalues and eigenvectors, are 
calculated. 

e LAPW2 (RHO): This stage derives valence densities from the obtained eigenvectors. 

e LCORE: It is responsible for computing the core states and their respective densities. 

e MIXER: In this step, a mixing process is performed for the input and output densities to enable 


further iterations. 


The diagram (Figure 1) clearly illustrates how the various WIEN2K programs are used and 


executed and provides a visual orientation for better understanding. 
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Figure.1 The flowchart shows the program flow in the WIEN2K package. 
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B.2 VASP Code 


2.1 Main Input/Output Files in VASP 


VASP searches the specified directory for four primary input files: POSCAR, INCAR, KPOINTS, 
and POTCAR. 


Within the POSCAR file we specify the placement of a single atom within a box. This file 


encapsulates the initial grid geometry and ion positions. 


The INCAR file contains tags that control the calculation. These tags are detailed in the VASP 
wiki under specific categories and have default values if not explicitly defined in the INCAR 
file. As a central input file, INCAR includes a significant portion of the key words that are 
essential for the calculations. This includes parameters such as the limit energy, smearing 


parameters, convergence criteria and more. 


The KPOINT file specifies the coordinates and associated weights of k-points within the 
Brillouin zone for sampling. It also contains information about how these k-points are created, 


including the method used, such as the Monkhorst-Pack scheme. 


The POTCAR file contains pseudopotential data and associated information for the designated 


atoms. 


Regarding the output files, VASP generates three main files: CONTCAR, OUTCAR and 
OSZICAR. 


e The CONTCAR includes the optimized grid geometry and ion positions. 


e The OUTCAR file serves as the central output document in VASP and includes a 


comprehensive set of data generated during the calculation. 


e The OSZICAR file provides an optimized overview of each electronic step, presenting 
important information on a single line. This includes the number of iterations, total energy 


and variation in total energy. 


195 


Appendices 


2.2 Running Calculations 


Similar to the WIEN2K code, the VASP code is based on two central phases: cell relaxation and 


the SCF cycle. When run with the vasp_std executable, VASP can be run in parallel with mpirun. First, 


the electronic charge density is determined from the POTCAR file and remains constant for the first 


steps. This is attributed to the iterative diagonalization of the Hamiltonian, meaning that the Kohn-Sham 


(KS) orbitals used to update the charge density require a warm-up period to be considered reliable. 


Figure 2 provides another visualization of the SCF cycle and OUTPUT file generation. 
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Figure.2 Schematic representation of the self-consistent loop for solving the Kohn-Sham equation. 


(For the two spins, it is necessary to go through two such loops at the same time.) 
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B.3 Possibilities with WIEN2K and VASP Codes 


The WIEN2K and VASP codes are versatile tools for studying a wide range of properties in periodic 
materials. They allow researchers to determine important structural attributes such as structure type, 
lattice constants, bonding, bond angles, bulk modulus, formation and cohesion energies, phonon 
spectrum and molecular dynamics. They also enable in-depth study of electronic properties, including 
spin-resolved electronic band structures, densities of states (both total and projected) and magnetic 
moments. In addition, they facilitate the assessment of real and imaginary components of dielectric 


function as well as the evaluation of mechanical, thermal and electronic transport properties. 


e Structural Properties 


Ensuring a well-defined and relaxed geometry is a crucial aspect when examining crystals or 
surfaces. The aim is to find the geometry with the lowest energy. Since we are studying periodic systems, 
two very important concepts come into play: the crystal unit cells and the positions of ions, both of 
which define the crystal structure. Within WIEN2K, optimization is carried out using the eosfit program. 
This tool allows varying grid parameters up to a certain percentage. One can then graph the energy 
versus volume and fit this nonlinear relationship using various equations of state. WIEN2K provides 
options such as the Murnaghan [8] and Birch Murnaghan [9] equations of state. VASP offers various 
relaxation methods. Among these, the most commonly used approach is to relax the ion positions while 
keeping the cell shape constant, a configuration achievable with ISIF=4. Visualization was facilitated 


by Xcrysden [10] for WIEN2K and Vesta [11] for VASP code. 


To deal with phonon spectra, an additional package called PHONOPY (included in VASP) was 
used [12]. This utility constructs a supercell to calculate force constants in real space using DFPT 
(Density Functional Perturbation Theory) [13]. The PHONOPY code is used to calculate the phonon 


spectrum, a process performed for all monolayers studied in this dissertation. 


By applying ab-initio molecular dynamics simulations (AIMD), a better assessment of thermal 
stability can be achieved. AIMD is a computational approach that combines quantum mechanical 
principles with classical molecular dynamics to simulate the dynamic behavior of atoms and molecules 
in real time. Without going further into the mathematical equations that describe their method, let us 
simply note that the thermostat principle is to adjust the velocities of the system particles during the 
simulation so that the average kinetic energy of the system is equal to the kinetic energy corresponding 
to the target temperature. Typical AIMD simulations are limited to systems with a few hundred atoms 
and a total duration of ten to a thousand picoseconds [14]. To run a basic AIMD simulation in VASP 


(as in our case for monolayers), a simple copy of the input files along with changes to a few flags and 
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selecting the correct thermostat type is enough. It should be noted that the choice of thermostats to use 
depends on the specific requirements of the simulation and the accuracy required for the study being 


carried out. Various thermostats can be used to control temperature, including: 
1. Nose-Hoover Thermostat : 


This is a popular choice for maintaining a constant temperature during AIMD simulations. It uses a 


number of additional degrees of freedom to control the kinetic energy of the system. 
2. Langevin Thermostat : 


The Langevin thermostat adds a stochastic term to the equations of motion, simulating the effects of a 


heat bath. This allows the system to reach and maintain a desired temperature. 
3. Andersen Thermostat : 


The Andersen thermostat randomly selects particles and gives them a velocity that corresponds to a 


Maxwell-Boltzmann distribution at the desired temperature. 


In this dissertation, we implemented the Nose-Hoover thermostat. 


e Electronic Properties 


After geometry relaxation, the electronic density of states (DOS) was determined by integrating the 
electron density over k-space. Another important electronic property, the band structure, was also 
estimated by locating the high symmetry points in the Brillouin zone of the solid. In our WIEN2K 
analysis, we used the interface tailored for the post-processing of first principles data, equipped with the 
essential programs for property calculation. Meanwhile, in our work with VASP, we used the VaspKit 
package [15] for data post-processing. VaspKit proves to be an accessible toolkit that streamlines the 
initial setup of calculations and enables detailed post-processing analysis to derive a variety of material 
properties from the raw data generated by the VASP code. (All upcoming properties are determined 


using the same scheme). 


e Magnetic Properties 


Evaluating magnetic properties is crucial for applications such as spintronic devices and spin-based 
devices. These properties have been thoroughly investigated by analyzing spin-polarized density of 
states diagrams and spin-resolved electronic band structures, magnetic moments, and exchange 


constants. 
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e Mechanical Properties 


Mechanical properties are the properties of a material that describe how it reacts to applied forces 
or loads. These properties are important for understanding how a material behaves under different 
conditions, which is crucial for the design and construction of various structures and components. The 
most important mechanical properties include stiffness, hardness, ductility and brittleness, all of which 
can be determined from elastic constants. In Figure 3 we have shown the strain curve (g), which 
illustrates the deformation of the material in relation to the applied stress (6). Two main regions are 
identified: the first is called the “elastic region” and is characterized by a linear variation of deformation 
with stress, where the deformation is reversible; when stress is removed, the material returns to its 
original shape. The second region is the “plastic region” in which the deformation is irreversible. The 
stress value at which the plastic region begins is called the yield point. Note that the elastic and plastic 
regions correspond to low and high stresses, respectively. In this dissertation, we will work in the elastic 


range. 
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Figure.3 Schematic Stress-Strain Curve. After Ref [16]. 


In the linear elastic regime, the relationship between the stress (o) and the applied strain (¢) in solids 


follows the generalized Hooke's law [16] and can be expressed using Voigt notation [17]: 


Gi, = Ljat C iE; (eq.20) 
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The strain, stress, and elastic constants are each represented as tensor of 6x6 independent 
components denoted by subscripts i and j. Each code handles the elastic constants in its own way. 
WIEN2K typically begins with a perfect crystal lattice structure. To calculate elastic constants, the 
program applies small deformations (strains) along various directions to the crystal lattice. For each 
strain applied, WIEN2K calculates the total energy of the deformed crystal structure. By fitting these 
energy-strain curves, it is possible to extract the elastic constants and determine all elastic moduli. 
Similar to WIEN2K, VASP can calculate elastic constants by applying small deformations to the crystal 
lattice. This is often done using the second derivative of the total energy according to the atom positions, 
the so-called Hessian matrix. This matrix contains information about the curvature of the energy surface, 


which is related to the elastic constants. 


e Optical Properties 


It is crucial to understand how materials interact with light, and their optical properties provide a 
valuable tool for this. These properties are characterized by the complex dielectric function ¢(@), which 


determines the response of the medium to the electromagnetic field [18]: 


E(W) = €4(@) + i€2(@) (eq.21) 


Here, €, and €2 represent the real and imaginary parts of the complex dielectric function, respectively, 


q@ denotes the photon frequency. 


In the WIEN2K framework, the "optic" subroutine is employed to compute the optical properties. 
This subroutine calculates the components of the matrix dipole moment for each k-point and each 
combination of occupied and empty bands. The determination of €2 components across the Brillouin 
zone is accomplished by the "joint" subroutine, generating the file case.joint. Application of the 
Kramers-Kronig formula to compute the €, components is carried out by the "kram" subroutine. At 
this stage, we specify the value of the "scissor operator," determined by the difference between the 
experimentally measured optical gap and the theoretically derived one. The resulting files include 


case.epsilon, case.sigmak, case.absorption, case.reflection, case.refraction, and case.eloss. 


Within the VASP code, the VASPKIT package proves invaluable in determining frequency- 
dependent dielectric functions. This package facilitates the extraction of linear spectra by parsing the 
output files using FLAG 71. The resulting files are labeled with names corresponding to each specific 
property, including ABSORPTION_2D, OPTICAL CONDUCTIVITY _2D, REFLECION_2D, and 
TRANSMISSION_2D. 
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e Thermoelectric Properties 


Thermoelectric properties refer to the characteristics of a material that govern its ability to convert 
temperature gradients into electrical voltage or vice versa, a phenomenon known as the Seebeck effect. 
These properties are crucial in thermoelectric applications, which involve harnessing waste heat or 
providing localized cooling. Optimizing the Figure of merit (ZT) for thermoelectric performance is a 
key goal in thermoelectric materials research. This parameter depends on properties like electrical 
conductivity (6), thermal conductivity (K= Ke+ k1), and the Seebeck coefficient (S) that could be 
manipulated through simulation to enhance the overall performance. The thermoelectric properties are 


calculated using the BoltzTrap2 code [19]. 
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